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We investigate weakly coupled spin-1/2 ladders in a magnetic field. The work is motivated by 
recent experiments on the compound (C 5 Hi2N) 2 CuBr4 (BPCB). We use a combination of numerical 
and analytical methods, in particular the density matrix renormalization group (DMRG) technique, 
to explore the phase diagram and the excitation spectra of such a system. We give detailed results 
on the temperature dependence of the magnetization and the specific heat, and the magnetic field 
dependence of the nuclear magnetic resonance (NMR) relaxation rate of single ladders. For cou- 
pled ladders, treating the weak inter ladder coupling within a mean- field or quantum Monte Carlo 
approach, we compute the transition temperature of triplet condensation and its corresponding an- 
tiferromagnetic order parameter. Existing experimental measurements are discussed and compared 
to our theoretical results. Furthermore we compute, using time dependent DMRG, the dynamical 
correlations of a single spin ladder. Our results allow to directly describe the inelastic neutron scat- 
tering cross section up to high energies. We focus on the evolution of the spectra with the magnetic 
field and compare their behavior for different couplings. The characteristic features of the spectra 
are interpreted using different analytical approaches such as the mapping onto a spin chain, a Lut- 
tinger liquid (LL) or onto a t-J model. For values of parameters for which such measurements exist, 
we compare our results to inelastic neutron scattering experiments on the compound BPCB and 
find excellent agreement. We make additional predictions for the high energy part of the spectrum 
that are potentially testable in future experiments. 

PACS numbers: 75.10.Jm, 75.40.Gb, 75.40.Cx, 75.30.Kz 



I. INTRODUCTION 



Many fascinating magnetic properties of solids are re- 
lated to quantum effect si. In particular, due to the Pauli 
principle, the interplay between interactions and kinetic 
energy can induce a strong antiferromagnetic spin ex- 
change. Such exchange leads to a remarkable dynamics 
for the spin degrees of freedom. On simple structures, 
the antiferromagnetic exchange can stabilize an antifer- 
romagnetic order. By variations in dimensionality and 
connectivity of the lattice a variety of complex phenom- 
ena can arise. 

Recently, among those effects two fascinating situa- 
tions in which the interaction strongly favors the forma- 
tion of dimers have been explored in detail. The first 



situation concerns a high dimensional system in which 
the antiferromagnetic coupling can lead to a spin liq- 
uid state made of singlets along the dimers. In such a 
spin liquid the application of a magnetic field leads to 
the creation of triplons which are spin-1 excitations. The 
triplons which behave essentially like itinerant bosons can 
condense leading to a quantum phase transition that is 
in the universality class of Bose-Einstein condensation^— 
(BEC). Such transitions have been explored experimen- 
tally and theoretically in a large variety of materials, be- 
longing to different structures and dimensionalities^. On 
the other hand, low dimensional systems behave quite 
differently. Quantum fluctuations are extreme, and no 
ordered state is usually possible. In many quasi one- 
dimensional systems the ground state properties are de- 
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scribed by Luttinger liquid (LL) physics 6,7 that predicts 
a quasi long range order. The elementary excitations are 
spin- 1/2 excitations (spinons). They behave essentially 
as interacting spinless fermions. This typical behaviour 
can be observed in spin ladder systems in the presence 
of a magnetic field. Although such systems have been 
studied theoretically intensively for many years in both 
zerc£— and finite magnetic field 2 -^ - — , a quantitive de- 
scription of the LL low energy physics remained to be 
performed specially for a direct comparison with experi- 
ments. 

Quite recently the remarkable ladder compound 28 
(C 5 Hi 2 N) 2 CuBr 4 , usually called BPCB (also known as 
(Hpip)2CuBr4), has been investigated. The compound 
BPCB has been identified to be a very good realiza- 
tion of weakly coupled spin ladders. The fact that the 
interladder coupling is much smaller than the intralad- 
der coupling leads to a clear separation of energy scales. 
Due to this separation the compound offers the exciting 
possibility to study both the phase with Luttinger liq- 
uid properties typical for low dimensional systems and 
the BEC condensed phase typical for high dimensions. 
Additionally, the magnetic field required for the realiza- 
tion of different phases lies for this compound in the ex- 
perimentally reachable range. The LL predictions have 
been quantitatively tested for magnetization and specific 
heati22, nuclear magnetic reasonance 3 -^ (NMR) and neu- 
tron diffraction^ measurements. Additionally the BEC 
transition and its corresponding order parameter have 
experimentally been observed by NMR 30 and neutron 
diffraction measurements 3 -^. 

In addition, the excitations of this compound have re- 
cently been observed by inelastic neutron scatterin g 32 ! 33 
experiments (INS). These are directly related to the dy- 
namical correlations of spin ladders in a magnetic field. 
These dynamical correlations have hardly been investi- 
gated so far. The direct investigation of such excitations 
is of high interest, since they not only characterize well 
the spin system, but the properties of the triplon/spinon 
excitations are also closely related to the properties of 
some itinerant bosonic/fermionic systems. Indeed using 
such mappings 6 of spin systems to itinerant fermionic or 
bosonic systems, the quantum spin systems can be used 
as quantum simulators to address some of the issues of 
itinerant quantum systems. One of their advantage com- 
pared to regular itinerant systems is the fact that the 
Hamiltonian of a spin system is in general well charac- 
terized, since the spin exchange constants can be directly 
measured. The exchange between the spins would corre- 
spond to short range interactions, leading to very good 
realization of some of the models of itinerant particles, for 
which the short range of the interaction is usually only 
an approximation. In that respect quantum spin systems 
play a role similar to the one of cold atomic gases 3 -^, in 
connection with the question of itinerant interacting sys- 
tems. 

In this paper, we present a detailed calculation of the 
properties of weakly coupled spin- 1/2 ladders. We fo- 



cus in particular on their dynamics and their low energy 
physics providing a detailled analysis and a quantitative 
description necessary for an unbiased comparison with 
experiments. More precisely, we explore the phase dia- 
gram of such a system computing static quantities (mag- 
netization, specific heat, BEC critical temperature, order 
parameter) and the NMR relaxation rate, using a com- 
bination of analytic (mostly Luttinger liquid theory and 
Bethe ansatz (BA)) and numerical (mostly density ma- 
trix renormalization group (DMRG) and quantum Monte 
Carlo (QMC)) techniques. We compare our results with 
the various measurements on the compound BPCB. A 
short account of some of these results in connection 
with measurements on BPCB was previously published 
in Refs. l29l - [3ll We here extend these results and give the 
details on how the theoretical results were obtained. Mo- 
tivated by recent experimental measurements, we further 
investigate the excitation spectra and dynamical corre- 
lation functions at high and intermediate energies, for 
which a theoretical description is very challenging. We 
show how for the low energy part of the spectrum it is 
possible to use the mapping to low energy effective the- 
ories such as the LL or to a spin chain which can be 
solved by Bethe ansatz technique s 35 ! 36 . Such a technique 
does not work, however, for energies of the order of the 
magnetic exchange of the system. In this manuscript we 
thus complement such analytical approaches by a DMRG 
analysis. We use the recent real-time variant to obtain 
the dynamics 3 - 7 . — in real time and the dynamical corre- 
lation functions. The same technique can also be used 
to obtain finite temperature results^"—. This allows to 
obtain an accurate computation of the excitation spec- 
tra and correlation functions in the high energy regime. 
We use different analytical approaches to interpret our 
numerical results. 

The plan of the paper is as follows. Sec. [II] defines 
the model of weakly coupled spin ladders. Its basic ex- 
citations and phase diagram are introduced as well as 
the spin chain mapping which proves to be very helpful 
for the physical interpretations. Sec. IIIII briefly recalls 
the different analytical (LL, BA) and numerical (DMRG, 
QMC) techniques which we used to obtain the results 
described in the present paper. Sec. IIVI gives a detailed 
characterization of the phase diagram focusing on the 
static properties and the NMR relaxation rate. Sec. [V] 
presents the computed dynamical correlations of a sin- 
gle spin ladder at different magnetic fields and couplings. 
The numerical calculations are compared to previous re- 
sults (link cluster expansion, spin chain mapping, weak 
coupling approach) and analytical descriptions (LL, t-J 
model). Sec. IVII directly compares the computed quan- 
tities to experimental measurements. In particular the 
theoretical spectra are compared to the low energy INS 
measurements on the compound BPCB. It also provides 
predictions for the high energy part of the INS cross sec- 
tion. Finally, Sec. IVIII summarizes our conclusions and 
discusses further perspectives. 



II. COUPLED SPIN-1/2 LADDERS 

In this section we introduce the theoretical model of 
weakly coupled spin- 1/2 ladders in a magnetic field. We 
recall its low-temperature phase diagram, paying spe- 
cial attention to the regime of strong coupling along the 
rungs of the ladder. This regime is particularly inter- 
esting since it is realized in the spin-ladder compound 
(C 5 Hi2N) 2 CuBr4, customarily called BPCB. We discuss 
briefly the energy scales for BPCB in the present section, 
leaving more detailled discussions for Sec. EH 



A. Model 

The Hamiltonian we consider is 

Hsb = + J' ^ &l,k,n ' Si',k',ii f - (1) 

Here is the Hamiltonian of the single ladder \i and 
J' is the strength of the interladder coupling. The op- 
erator S lik ,p = (Si^Sf^Sf^) acts at the site I 
(I = 1,2,...,L) of the leg k (k = 1,2) of the lad- 
der \i. Often we will omit ladder indices from the sub- 
scripts of the operators (in particular, replace with 
Si y k) to lighten notation. Sj* k (a = x,y,z) are conven- 
tional spin- 1/2 operators with [Sf k ,Sf k ] = iSf k , and 



^i,k 

The Hamiltonian of the spin- 1/2 two- leg ladder 
illustrated in Fig. [T] is 



Ha = J±H± + J\\H\\ 



(2) 



where J± (Jn) is the coupling constant along the rungs 
(legs) and 



i 

H\\ = ^2 ^ l ' k ' ^l+l,k 



(3) 
(4) 



The magnetic field, h z , is applied in the z direction, 
and M z is the z-component of the total spin opera- 
tor M = + S; } 2). Since has the symmetry 
h z -h z , M z U —M z , we only consider h z > 0. The 
relation between h z and the physical magnetic field in 
experimental units is given in Eq. ([40]) . 



B. Energy scales 

In the present paper we focus on the case of spin- 1/2 
antiferromagnetic ladders weakly coupled to one another. 
This means that the interladder coupling J' > is much 
smaller than the intraladder couplings Jm and J_l, i.e. 



FIG. 1. Single ladder structure: J± (Jy) is the coupling along 
the rungs (legs) represented by thick (thin) dashed lines and 
Si : k are the spin operators acting on the site / of the leg 
k = 1,2. 



As we will show, the model (pp) accurately describes the 
magnetic properties of the compound BPCB. Its detailed 
description is given in Sec. EH The couplings have been 
experimentally determined to b o 3Q i 31 



and 30 



J' « 20 - 100 mK 



J\\ « 3.55 K, J ± « 12.6 K. 



(6) 



(7) 



More details about the determination of the couplings for 
the compound BPCB are given in Sec. IVT1 



C. Spin ladder to spin chain mapping 

The physical properties of a single ladder ((2]) are de- 
fined by the value of the dimensionless coupling 



J\\ 

1 = T,- 



(8) 



In the limit Jm = (therefore 7 = 0) the rungs of the lad- 
der are decoupled. We denote this decoupled bond limit 
hereafter. The four eigenstates of each decoupled rung 
are: the singlet state 



\s) 



\n) - lit) 

V2 



(9) 



with the energy E s = — 3J^/4, spin 5 = 0, and z- 
projection of the spin S z = 0, and three triplet states 



Itt), \t°) 



in) + lit) 

V2 ' 



|t-> = \U) (io) 



< J' < Jm and J_l. 



(5) 



with S = 1, S z = 1,0, —1, and energies E t + = J±/4 — h z , 
E t o = J_l/4, E t - = J_l/4 + h z , respectively. The ground 
state is \s) below the critical value of the magnetic field, 
h^ BL = J_l, and above. The dependence of the 
energies on the magnetic field is shown in Fig. [2ja. 

A small but finite 7 > delocalizes triplets and creates 
bands of excitations with a bandwidth ~ Jn for each 
triplet branch. This leads to three distinct phases in the 
ladder system (|2]) depending on the magnetic field: 

(i) Spin liquid phased, which is characterized by a 
spin-singlet ground state (see Sec. IIV A|) and a 
gapped excitation spectrum (see Sec. IV Bp . This 
phase appears for magnetic fields ranging from 
to h c , . 




liquid / \ LL / \ polarized 



/ Quantum \ 
critical \ 

T (b) 

FIG. 2. (a) Energy of the triplets |t + ), |t°), \t~) (solid lines) 
and singlet \s) (dashed line) versus the applied magnetic field 
in the absence of an interrung coupling (Jy = 0). The dot- 
ted lines represent the limits of the triplets excitation band 
when J|| ^ 0. (b) Phase diagram of weakly coupled spin lad- 
ders: crossovers (dotted lines) and phase transition (solid line) 
that only exists in the presence of an interladder coupling are 
sketched. 



with the new effective spin- 1/2 operators Sf. 

S tk = ^f^' Sl k = \(l + 2~Sf). (12) 

The Hamiltonian (|2j) reduces to the Hamiltonian of the 
spin- 1/2 XXZ Heisenberg chain 

fez = J\\ (Sf 5f+i + SfSf +1 + ASi Sf +1 ) 

I 

-^^(-T + f-y)- < 13 > 

Here the pseudo spin magnetization is M z = ^ Sf, the 
magnetic field h z = h z — J± — Jy/2 and the anisotropy 
parameter 



Note that the spin chain mapping constitutes a part of a 
more general strong coupling expansion of the model (|2]), 
as discussed in the appendix [A] 

For the compound BPCB the parameter 7 is rather 
small 

7 « 3^ « 0.282. (15) 

and the spin chain mapping (fT3|) gives the values of many 
observables decently well. Some important effects are, 
however, not captured by this approximation. Examples 
will be given in later sections. 



(ii) Gapless phase, which is characterized by a gapless 
excitation spectrum. It occurs between the critical 
fields h Cl and h C2 . The ground state magnetization 
per rung, m z = (M z )/L, increases from to 1 for 
h z running from h Cl to h C2 . The low energy physics 
can be described by the LL theory (see Sec. IIII C|) . 

(iii) Fully polarized phase, which is characterized by the 
fully polarized ground state and a gapped excita- 
tion spectrum. This phase appears above h C 2- 

Besides ladders the transition between (i) and (ii) can 
occur in several other gapped systems such as Haldane 
S = 1 chains or frustrated chai mj 19 i 45 ~ — . In the gapless 
phase, the distance between the ground state and the 
bands \t°) and \t~), which is of the order of J_l, is much 
larger than the width of the band ~ Jm , since 7 <C 1. 

For small 7 the ladder problem can be reduced to a 
simpler spin chain problem. The essence of the spin chain 
mappingl i 17 i 48 i 49 is to project out \t°) and \t~) bands 
from the Hilbert space of the model (|2]). The remaining 
states I s) and are identified with the spin states 

ll> = l*>, lf> = l* + >- (11) 

The local spin operators can therefore be identified 
in the reduced Hilbert space spanned by the states JTT} 



D. Role of weak interladder coupling 

Let us now turn back to the more general Hamilto- 
nian (pQ) and discuss the role of a weak interladder cou- 
pling J' (couplings ordered as in Eq. fl5J). The spin liq- 
uid and fully polarized phases are almost unaffected by 
the presence of J' whenever the gap in the excitation 
spectrum is larger than J' (see, e.g., Ref. fl8l for more de- 
tails). However, a new 3D antiferromagnetic order in the 
plane perpendicular to h z emerges in the gapless phase 
for T < J' . The corresponding phase, called 3D- ordered, 
shows up at low enough temperatures T c in numerous 
experimental systems with reduced dimensionality and a 
gapless spectrum 5 . For the temperature T > J' the lad- 
ders decouple from each other and the system undergoes 
a deconfinement transition into a Luttinger liquid regime 
(which will be described in Sec. HITC]) . For T > J\\ the 
rungs decouple from each other and the system becomes 
a (quantum disordered) paramagnet. All the above men- 
tioned phases are illustrated in Fig. [2jb. 

III. METHODS 

In this section, we present the methods used to study 
the ladder system and its mean-field extension to the 
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case of weakly coupled ladders. We first focus on the so 
called density matrix renormalization group (DMRG) or 
matrix product state (MPS) methods. These numerical 
methods allow us to investigate dynamical correlations 
at zero and finite temperature. Additionally we discuss 
the Bethe ansatz used to obtain properties of the system 
after the spin chain mapping. Furthermore we introduce 
an analytical low energy description for the gapless phase, 
the Luttinger liquid theory. This theory in combination 
with a numerical determination of its parameters (see 
appendix [Bj gives a quantitative description of the low 
energy physics. Finally, we treat the weak inter ladder 
coupling J' by a mean field approach, both analytically 
and numerically, and a quantum Monte Carlo (QMC) 
technique. 



A. DMRG 

A numerical method used to determine static and dy- 
namical quantities at zero and finite temperature of a 
quasi one-dimensional system is the DMRG. This method 
was originally introduced by S.R. White 50 to study static 
properties of one dimensional systems. Since usually the 
dimension of the total Hilbert space of a many-body 
quantum system is too large to be treated exactly, the 
main idea of the DMRG algorithm is to describe the im- 
portant physics using a reduced effective space. This re- 
duced effective space is chosen optimally by using a vari- 
ational principle. The DMRG has been proven very suc- 
cessfully in many situations and has been generalized to 
compute dynamical properties of quantum systems using 
different approaches in frequency spaced— . Recently 
the interest in this method even increased after a suc- 
cessful generalization to time-dependent phenomena and 
finite temperature situations^—. The real-time calcu- 
lations give an alternative route to determine dynamical 
properties of the system^ which we use in the follow- 
ing. An overview of the method, its extensions and its 
successful applications to real-time and finite tempera- 
ture can be found in Refs. [5l| and [HI. Further details 
on the method and its technical aspects are given in the 
appendix O 



the thermodynamic properties of the spin- 1/2 XXZ chain 
model ([T3|) by the Bethe ansatz technique is a simpler, 
but still non-trivial task, requiring a solution of an infi- 
nite set of non-linear coupled integral equations 56 . The 
solution of such equations can be only found numeri- 
cally, and already in the 1970s this was done with a high 
precision ' . 

Later on, an alternative to the Bethe ansatz, the quan- 
tum transfer matrix method, was used to get the thermo- 
dynamics of the XXZ chain in a magnetic fielc—. Within 
this approach the free energy of the system is expressed 
through the largest eigenvalue of the transfer matrix. 
This largest eigenvalue is given by the solution of a set 
of non- linear equations (Eq. (66. a) of Ref. l58l ) which are 
written in a form very suitable for solving them itera- 
tively. In the present paper we followed this route and 
got the results for the specific heat for all temperatures 
and various magnetic fields with very high precision, see 
Sec. lIVB21 



C. Luttinger Liquid (LL) 

The Luttinger liquid Hamiltonian governs the dynam- 
ics of the free bosonic excitations with linear spectrum 
and can be written as&i 



2tt/ 



dx 



uK (d x 0{x)Y 



(16) 



where <p and are canonically commuting bosonic fields, 
[(p(x),d y 0(y)] = iir5(x — y). The dimensionless param- 
eter K entering Eq. ([T6]) is customarily called the Lut- 
tinger parameter, and u is the propagation velocity of the 
bosonic excitations (velocity of sound). Many ID inter- 
acting quantum systems belong to the Luttinger Liquid 
(LL) universality class: the dynamics of their low-energy 
excitations is governed by the Hamiltonian ([T6]) and the 
local operators are written through the free boson fields 
and (j) (the latter procedure if often called bosonization). 

The spin- 1/2 XXZ Heisenberg chain, Eq. (fT3j) . in the 
gapless phase is a well-known example of a model belong- 
ing to the LL universality class. Its local operators are 
expressed through the boson fields as follows^: 



B. Bethe ansatz 

The spin- 1/2 XXZ chain ([T3|) which is obtained after 
the spin chain mapping of the system is exactly solvable: 
the so-called Bethe ansatz technique gives explicit ana- 
lytic expressions for its eigenfunctions and spectru m 54 ! 55 . 
To convert this information into a practical recipe of cal- 
culation of the correlation functions is a highly sophisti- 
cated problem. However, a known solution to this prob- 
lem (Ref. [35| and references therein) incorporates involved 
analytics and numerics, the latter limiting the precision 
of the final results to about the same extent as to-date 
implementations of the DMRG method. Calculation of 



S^x) 



\f2A x (-l) x 

+2^/B^cos(20(x) - 2nm z x) 



and 

S z (x) = m z 



d x (j){x) 



(17) 



y / 2A~ z (-l) x cos(20(x) - 27rm z x). (18) 



Here the continuous coordinate x = la is given in units of 
the lattice spacing a, rh z = (M z )/L is the magnetization 
per site of the spin chain, and A x , B x and A z are coeffi- 
cients which depend on the parameters of the model ([T3]) . 
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How to calculate K, u, A x , B x , and A z is described in 
the appendix IB 1[ 

The Hamiltonian (p~3|) is the leading term in the strong 
coupling expansion (the parameter 7 = J\\/J± <C 1, see 
appendix [A]) of the model (j2]). Local operators of the 
latter model are bosonized by combining Eqs. ([T2]h ([T7]h 
and JT5J). The analysis of the model ([2]) suggests ^ 19 ! 20 
that the bosonization of the local spins can be performed 
for any values of J± and J\\ in the gapless regime. We 
would like to stress that even for a small 7 some parame- 
ters out of K, u, A x , B X1 and A z show significant numer- 
ical differences if calculated within the spin chain ([T3]) 
compared to the spin ladder (|2j). We discuss this issue in 
the appendix IB II 



D. Mean- field approximation 

Up to now, we have presented methods adapted to 
deal with one dimensional systems. In real compounds, 
an interladder coupling is typically present. As discussed 
in Sec. IIID| in the incommensurate regime this inter- 
ladder coupling J' (cf. Eq. (pQ)) can lead to a new three 
dimensional order (3D-ordered phase in Fig. [2jb) at tem- 
peratures of the order of the coupling J'. In the case of 
BPCB the interladder coupling is much smaller than the 
coupling inside the ladders, i.e. / < Ji, Jn (Sec. IVI A]) . 
Therefore, unless one is extremely close to h& or h C 2 
one can treat the interladder coupling with a standard 
mean-field approximation. This approach incorporates 
all the fluctuations inside a ladder. However, it overes- 
timates the effect of J' by neglecting quantum fluctua- 
tions between different ladders. Such effects can partly 
be taken into account by a suitable change of the in- 
terladder coupling 31 that will be discussed in Sec. IIVDI 
Close to the critical fields the interladder coupling J' be- 
comes larger than the effective energy of the one dimen- 
sional system. This forces one to consider a three dimen- 
sional approach from the start and brings the physics 
of the system in the universality class of Bose-Einstein 
condensation 2 -^. In the following we consider that we are 
far enough (i.e. by an energy of the order of J') away 
from the critical points so that we can use the mean-field 
approximation. 

Since the single ladder correlation functions along the 
magnetic field direction (z-axis) decay faster than the 
staggered part of the ones in the perpendicular x?/-plane 
(see Eqs. (fB2|) and (lB3l) for the LL exponent K of the 
ladder shown in Fig. [23]), the three dimensional order will 
first occur in this plane. Thus the dominant order param- 
eter is the q = it staggered magnetization perpendicular 
to the applied magnetic field. The mean-field decoupling 
of the spin operators of neighboring ladders thus reads 



Sfk ~ -(-±) lJrk m x 



St k = \ ~ (-1) 



l+k m z 



(-l) 1+fc (^ fe ) (19) 
(20) 



We have chosen the ^-ordering to be along the x-axis 
and m z a will be very small and therefore neglected. 

This approximation applied on the interladder inter- 
action part of the 3D Hamiltonian Hsb (Eq. ([I])) leads 
to 



Hmf = J\\ H\\ 
n r J'm z 



,J'm x 



-£(-!) 

l,k 



l-\-k qx 



Sf k - (21) 



Here we assume that the coupling is dominated by n c 
neighboring ladders, where n c is the rung connectivity 
(n c = 4 for the case of BPCB, cf. Fig. [18]). This mean- 
field Hamiltonian corresponds to a single ladder in a 
site dependent magnetic field with a uniform component 
in the z-direction and a staggered component in the in- 
direction. The ground state wave function of the Hamil- 
tonian must be determined fulfilling the self-consistency 
condition for m z and m x using numerical or analytical 
methods. This amounts to minimize the ground state 
energy of some variational Hamiltonian. 



1. Numerical mean-field 

The order parameters m z and m x a can be computed 
numerically by treating the mean-field Hamiltonian Hmf 
self consistently with DMRG. These parameters are eval- 
uated recursively in the middle of the ladder (to min- 
imize the boundary effects) starting with m z = and 
m x = 0.5. An accuracy of < 10 -3 on these quantities 
is quickly reached after a few recursive iterations (typ- 
ically - 5) of the DMRG keeping few hundred DMRG 
states and treating a system of length L = 150. We veri- 
fied by keeping as well the alternating part of the z-order 
parameter ml that this term is negligible (< 10 -5 ). 



2. Analytical Mean-field 

Using the low energy LL description of our ladder sys- 
tem (see Sec. IIII C]h it is possible to treat the mean- field 
Hamiltonian Hmf within the bosonization technique. In- 
troducing the LL operators (fT7|) and (fl8|) in Hmf (f2T|) 
and keeping only the most relevant terms leads to the 
Hamiltonian^^ 



Hsg 



2tt/ 



dx 



uK (d x 0(x)f + ^ {d x (j>{x)) 2 



A x n c J'm x a / dxcos(0(x)) (22) 



where we neglected the mean-field renormalization of h z 
in (|2T|) . This Hamiltonian differs from the standard LL 
Hamiltonian J^ll ([16]) by a cosine term corresponding 
to the x-staggered magnetic field in ([2T]) . It is known 
as the sine- Gordon Hamiltonian ^ 61 1 62 . The expectation 
values of the fields can be derived from integrability 6 - 3 -. 
In particular m x can be determined self-consistently. 
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E. Quantum Monte Carlo (QMC) 

In order to take into account the detailled coupling 
structure of the BPCB compound shown in Fig. [TBI which 
is neglected in the mean-field approximation (Sec. IMP]) , 
we employ a stochastic series expansion implementation 
of the QMC technique with directed loop updates^ pro- 
vided with the ALPS librarie a 65 ^ 66 . Nevertheless due to 
the strong anisotropy of the couplings (j5j), the tempara- 
tures at which the effects of the interladder coupling J' 
become visible are not reachable with this method. The 
QMC results for the transition tem per ature of the 3D- 
ordered phase, T c , presented in Ref. |3jJ, Sec. II VP II and 
appendix [D] are then computed with a J' of ~ 3 times 
bigger than that extracted in Ref. [30| and Sec. IIVD II 
making the 3D effects numerically accessible. 



IV. STATIC PROPERTIES AND NMR 
RELAXATION RATE 

We begin our analysis of the different phases of the 
coupled spin ladder system, Fig. [2jb, by computing ther- 
modynamic quantities, such as the magnetization, the 
rung state density and the specific heat. In particular, we 
test the LL low energy prediction of the latter and evalu- 
ate the related crossover to the quantum critical regime. 
Furthermore we discuss the effect of the 3D interladder 
coupling computing the staggered magnetization in the 
3D-ordered phase and its critical temperature. We finally 
discuss the NMR relaxation rate in the gapless regime re- 
lated to the low energy dynamics. In order to compare 
these physical quantities to the experiments, all of them 
are computed for the BPCB parameters (see Sec. I VI A]) . 



A. Critical fields 

The zero temperature magnetization contains ex- 
tremely useful information. Its behavior directly gives 
the critical values of the magnetic fields h c \ and h C 2 at 
which the system enters and leaves the gapless regime, 
respectively (Fig. [2jb). In Fig. [3] the dependence of the 
magnetization on the applied magnetic field is shown for 
a single ladder and for the weakly coupled ladders. At low 
magnetic field, h z < h c \, the system is in the gapped spin 
liquid regime with zero magnetization, and spin singlets 
on the rungs dominate the behavior of the system^, see 
Fig. |4j At h z = h c i, the Zeemann interaction closes the 
spin gap to the rung triplet band (Fig. [2]). Above 
h z > h c i the triplet band starts to be populated 
leading to an increase of the magnetization with h z . The 
lower critical field in a 13th order expansion^ in 7 is 
h c i ~ 6.73 T for the BPCB parameters. At the same 
time the singlet and the high energy triplet occupation 
decreases, Fig. ffl For h z > h c2 = J± + 2Jy « 13.79 T 
(for the compound BPCB), the band is completely 



filled and the other bands are depopulated. The sys- 
tem becomes fully polarized (m z = 1) and gapped. The 
two critical fields, h c \ and h C 2, are closely related to the 
two ladder exchange couplings, J± and Jy. As they are 
experimentally easily accessible, assuming that a ladder 
Hamiltonian is an accurate description of the system, 
these critical fields can be used to determine the ladder 
couplings 30 . 

Such a general behavior of the magnetization is seen 
for both the single ladder and the weakly coupled lad- 
ders in Fig. [3j In particular, the effect of a small cou- 
pling J' between the ladders is completely negligible in 
the central part of the curve. Only in the vicinity of 
the critical fields, the single ladder and the coupled lad- 
ders show a distinct behavior. The single ladder be- 
haves like an empty (filled) one-dimensional system of 
non-interacting fermions which leads to a square-root be- 
havior m z oc (h z — h c \) x l 2 close to the lower critical field 
and 1 — m z oc (h C 2 — h z ) 1 / 2 close to the upper criti- 
cal field. In contrast, in the system of weakly coupled 
ladders, a 3D-ordered phase appears at low enough tem- 
peratures in the gapless regime (see Sec. Ill Dl and IIII D]) . 
The magnetization dependence close to the critical fields 
becomes linear, m z oc h z — h&, and 1 — m z oc h C 2 — h z , 
respectively 2 -^. In comparison with the single ladder, 
the critical fields given above are shifted by a value of 
the order of J 7 . This behavior is in the universality class 
of the Bose-Einstein condensation and is well reproduced 
by the mean-field approximation shown in the insets of 
Fig. [3] close to the critical fields. 

For comparison, the magnetization of a single lad- 
der in the spin chain mapping is also plotted in Fig. [3] 
This approximation reproduces well the general behavior 
of the ladder magnetization discussed above. However, 
note that for the exchange coupling constants considered 
here the lower critical field in this approximation is dif- 
ferent from the ladder one. The lower critical field is 
^i XZ = J± ~ J \\ ~ 6 - 34 T < fed. The upper critical 
field /^ xz = J_l + 2 J|| = h C 2 is the same as for the lad- 
der. If we rescale /i xxz and /^ xz to match the critical 
fields h cl and h c2 (h z -+ ^"fiT^l^ + h cl ), the 
magnetization curve gets very close to the one calculated 
for a ladder. However, in contrast to the magnetization 
curve for the ladder, the corresponding curve in the spin 
chain mapping is symmetric with respect to its center at 

7 XXZ I 7 xxz 

h xxz = n cl +n c2 = j ± + due to the absence of 
the high energy triplets. 



B. The Luttinger liquid regime and its crossover to 
the critical regime 

The thermodynamics of the spin- 1/2 ladders has been 
studied in the pas t n i 21 i 22 i 29 . We here summarize the 
main interesting features of the magnetization and the 
specific heat focussing on the crossover between the LL 
regime and the quantum critical region using the BPCB 
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FIG. 3. (Color online): Dependence of the magnetization per 
rung m z on the magnetic field h z at zero temperature for the 
single ladder with the BPCB couplings (Sec. IVI A[) (dashed 
red line), the spin chain mapping (dotted blue line) rescaled to 
fit with the single ladder critical fields (dash-dotted blue line), 
and for the weakly coupled ladders treated by the mean-field 
approximation (solid black line). The insets emphasize the 
different behavior of the magnetization curves for the single 
(dashed red line) and coupled (solid black line) ladders close 
to the critical fields which are indistiguishable in the main 
part of the figure. The dotted lines in the insets correspond 
to the linear and square root like critical behaviour. 
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FIG. 4. (Color online): Rung state density versus the applied 
magnetic field h z at zero temperature for the single ladder 
with the BPCB couplings (zero temperature DMRG calcu- 
lations averaging on the central sites of the ladder). The 
dash-dotted (black) lines correspond to the singlet density 
(p s )- The triplet densities are represented by the solid (red) 
lines for the dashed (blue) lines for (po) and the dotted 

(green) lines for {p~). 



parameters (Sec. IVI A]) . As the inter ladder exchange cou- 
pling J' is very small compared to the ladder exchange 
couplings J\\ and J_l, it is reasonable to neglect it in this 
regime far from the 3D phase. Therefore we focus on a 
single ladder in the following. 



1. Finite temperature magnetization 

We start the description of the temperature depen- 
dence of the magnetization in the two gapped regimes: 
the spin liquid phase and the fully polarized phase (not 
shown, cf. Refs. [2l| and |23[). For small magnetic fields 
h z < h c i, the magnetization vanishes exponentially at 
zero temperature and after a maximum at intermediate 
temperatures it decreases to zero for large temperatures. 
For large magnetic fields h z > h C 2 the magnetization in- 
creases exponentially up to m z = 1 at low temperature 
and decreases monotonously in the limit of infinite tem- 
perature. 

In the gapless regime, the magnetization at low tem- 
perature has a non-trivial behavior that strongly depends 
on the applied magnetic field. The temperature depen- 
dence of the magnetization computed with the T-DMRG 
(see appendix IC 2ft is shown in Fig. 0a for different values 
of the magnetic field in the gapless phase, h z = 9, 10, 11 T 
(h c i < h z < h c <i). In this regime new extrema appear in 
the magnetization at low temperature. This behavior 
can be understood close to the critical fields where the 
ladder can be described by a one-dimensional fermion 
model with negligible interaction between fermions. In- 
deed, in this simplified picture^ and in more refined 
calculation s 21 ! 22 ! 24 the magnetization has an extremum 
where the temperature reaches the chemical potential, 
i.e., at the temperature at which the energy of excita- 
tions starts to feel the curvature of the energy band. This 
specific behavior is illustrated in Fig. [5] a with the curve 
for h z = 11 T (h m = h <*+ h <* < h z < h c2 ). The low 
temperature maximum moves to higher temperature for 
h z < hm and goes over to the already discussed maximum 
for h z < h c \. Symmetrically with respect to h mi a low 
temperature minimum appears in the curve for h z = 9 T 
(h c i < h z < hm). This minimum slowly disappears for 
h z — )• h m (the curve for h z = 10 T is close to that). 

The location of the lowest extremum is a reasonable 
criterion to characterize the crossover temperature be- 
tween the LL and the quantum critical regim o 21 i 22 , since 
the extremum occurs at temperatures of the order of the 
chemical potential. A plot of this crossover temperature 
versus the magnetic field is presented in Fig. EJc. Follow- 
ing this criterium, the crossover has a continuous shape 
far from h m - Nevertheless, close to h m both extrema 
are close to each other (since the maximum still exists 
for h z < h m held at which the minimum appears). The 
criterium is thus not well defined. It presents a disconti- 
nuity at hm which is obviously an artefact. In the vicinity 
of /i m , we thus use another crossover criterium based on 
the specific heat (see Sec. IIVB 2|) that seems to give a 
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FIG. 5. (Color online): Temperature dependence of the 
magnetization per rung, m z (T), for (a) the ladder with the 
BPCB couplings (39]), and (b) the spin chain mapping at dif- 
ferent applied magnetic fields h z = 9 T (solid blue lines), 
h z = 10 T (dash-dotted green lines) and h z — 11 T (dashed 
red lines). The results were obtained using T-DMRG. The 
stars at T = K are the ground state magnetization per 
rung determined by zero temperature DMRG. The triangles 
(squares) marks the low (high) energy extrema. (c) Crossover 
temperature TLLof the LL to the quantum critical regime ver- 
sus the applied magnetic field (blue circles for the extremum 
in m z (T)\hz criterium and red crosses for the maximum in 
c(T)\h* criterium). 



more accurate description. In Ref. |29|, both criteria have 
been applied on the magnetocaloric effect and specific 
heat measurements on the compound BPCB, and give a 
crossover temperature in agreement with the computed 
ones. 

The temperature dependence of the magnetization of 
the spin chain mapping, Fig. Ob, exhibits a single low 
temperature maximum if h^ z < h z < h^ z (minimum 
if h^ z < h z < h^ z ). The appearance of a single ex- 
tremum and its convergence to m z — >> 0.5 when T — >> oo is 
due to the exact symmetry with respect to the magnetic 
field h^ z . This approximation reproduces the main low 
energy features of the ladder but fails to describe the high 
energy behavior which strongly depends on the high en- 
ergy triplets. 



2. Specific heat 

The specific heat has been investigated for similar pa- 
rameters as the ones considered in this paper in the 
gapped and gapless regime in Refs. [2l|, [23|, andl29l. Here 



we concentrate on the detailed analysis in the gapless 
regime, in particular on the low temperature behavior 
and the determination of a crossover temperature from 
the first maximum. We show in Fig. [6] the typical tem- 
perature dependence of the specific heat for different val- 
ues of the magnetic field. Comparisons with actual ex- 
perimental dataS for BPCB are excellent (see Fig. [6jb) . 
For these comparisons the theoretical data are computed 
with g = 2.06 related to the experimental orientation 
of the sample with respect to the magnetic fields (see 
Sec. I VI A|) and rescaled by a factor 0.98 in agreement 
with the global experimental uncertainties^. 

At low temperatures the specific heat has a contribu- 
tion due to the gapless spinon excitations. This results in 
a peak around 1.5 K. This peak is most pronounced 
for the magnetic field values lying mid value between the 
two critical fields. At higher temperatures the contribu- 
tion of the gapped triplet excitations leads to a second 
peak whose position depends on the magnetic field. To 
separate out the contribution from the low lying spinon 
excitations, we compare the specific heat of the ladder to 
the results obtained by the spin chain mapping in which 
we just keep the lowest two modes of the ladder (see 
Sec. Ill CI and appendix [A]) . The resulting effective chain 
model is solved using Bethe-ansatz and T-DMRG meth- 
ods. The agreement between these methods is excellent 
and the corresponding curves in Fig. [6] can hardly be dis- 
tinguished. However, a clear difference with the full spin 
ladder results is revealed. While at low temperatures the 
curves are very close, the first peak in the spin chain map- 
ping already lacks some weight, which stems from higher 
modes of the ladder. 

In Fig. [7] the low temperature region is analyzed in 
more detail. At very low temperatures the spinon modes 
of the ladder can be described by the LL theory (see 
Sec. IIII C|) which predicts a linear rise with temperature 
inversely proportional to the spinon velocity^^ 



cll(T) = 



Ttt 
3u ' 



(23) 



In Fig. [71 we compare the results of the LL, the Bethe 
ansatz and the DMRG results for the effective spin chain 
and the numerical DMRG results taking the full ladder 
into account. The numerical results for the adaptive T- 
DMRG at finite temperature are extrapolated to zero 
temperature by connecting algebraically to zero temper- 
ature DMRG results (see appendix IC 2|) . A very good 
agreement between ([23]) and numerics is found for low 
temperatures. However, at higher temperatures, the 
slope of the T LL description slightly changes with 
respect to the curves calculated with other methods. 
This change of slope reflects the fact that the curvature 
of the energy dispersion must be taken into account when 
computing the finite temperature specific heat, and this 
even when the temperature is quite small compared to 
the effective energy bandwidth of the system. The effec- 
tive spin chain and the numerical results for the ladder 
agree for higher temperatures (depending on the mag- 
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FIG. 6. (Color online): Specific heat per rung c versus the 
temperature T for different applied magnetic fields in the gap- 
less regime, (a) Full ladder T-DMRG calculations with the 
BPCB couplings J39} for h z = 9 T (solid blue line), h z = 10 T 
(dash-dotted green line), and h z = 11 T (dashed red line). 
Spin chain mapping at h z = 10 T solved by DMRG (dashed 
black line) and by Bethe ansatz (dotted black line). Note 
that the two lines are hardly distinguishable. The triangles 
(squares) mark the low (high) energy maxima of the specific 
heat versus temperature. The vertical dashed line marks the 
temperature T — 0.4 K below which the DMRG results are ex- 
trapolated (see appendix IC 2[) . (b) Comparison between mea- 
surements on the compound BPCB from Ref. 29 (dots) and 
the T-DMRG calculations^ (solid lines) at (b.l) h z = 9 T, 
(b.2) h z = 10 T and (b.3) h z = 11 T. 



netic field), before the higher modes of the ladder cause 
deviations. 

As for the magnetization (Sec. IIVB l|h the location 
of the low temperature peak can be interpreted as the 
crossover of the LL to the quantum critical regime. In- 
deed, in a free fermion description which is accurate 
close to the critical fields, this peak appears at the tem- 
perature for which the excitations stem from the bot- 
tom of the energy band. The corresponding temperature 
crossover is compared in Fig. 0c to the crossover temper- 
ature extracted from the first magnetization extremum 
(Sec. IIVB 1[) . The two crossover criteria are complemen- 
tary due to their domain of validity. The first specific 
heat maximum is well pronounced only in the center of 
the gapless phase. In contrast in this regime the pres- 
ence of two extrema close to each other in the magnetiza- 
tion renders the magnetization criterium very imprecise 
(cf. Sec. IIVB ID . 



C. Spin-lattice relaxation rate 
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FIG. 7. (Color online): Low temperature dependence of the 
specific heat per rung versus the temperature, c(T), for (a) 
h z = 9 T, (b) h z = 10 T, (c) h z = 11 T. The T-DMRG cal- 
culations are in red thick lines for the ladder with the BPCB 
couplings (|39p (black thick lines for the spin chain mapping). 
The two curves can hardly be distinguished. Their low tem- 
perature polynomial extrapolation is plotted in thin lines be- 
low T = 0.4 K (represented by a vertical dashed line). The 
linear low temperature behavior of the LL is represented by 
dashed lines (red for the ladder, black for the spin chain map- 
ping) . The dashed yellow lines correspond to the Bethe ansatz 
computation for the spin chain mapping. 



function! x * x ( x = 0,t) defined as in Eq. (tB6l) 



Im 



UJ 



dte lUJot X x a x (x = 0,t) (24) 



where T ^> ujq is assumed with the Larmor frequency ujq. 
7 n = 19.3 MHz/T is the nuclear gyromagnetic ratio of 
the measured nucleus ( 14 N in Ref. |3(J for BPCB) and A± 
is the transverse hyperfine coupling constant. 

Assuming Ju ^> T, the T x _1 due to the spin dynamics 
can be computed in the gapless regime using the LL low 
energy description. Following Ref. [2| we introduce the LL 
correlation (|B7|) into Eq. (fM|h and obtain 



Tf 1 



r 



l AjA x cos(^) (2kT 



B 



1 

4K' 



1 

2K 
(25) 

According to the known LL parameters (Fig. [23]) the 
shape of T{ l (h z ) plotted in Fig. [8] at T = 250 mK > T c 
is strongly asymmetric with respect to the middle of 
the gapless phase. The only free (scaling) parameter, 
A± = 0.057 T, is deduced from the fit of Eq. fl25|) to 
the experimental dataSi and agrees with other measure- 
ments. For comparison, the T x _1 obtained in the spin 
chain mapping approximation is also plotted in Fig. [8] 
As for other physical quantities, this description fails to 
reproduce the non-symmetric shape. 



As the spin-lattice relaxation in quantum spin sys- 
tems is due to pure magnetic coupling between electronic 
and nuclear spins, the NMR spin-lattice relaxation rate 
is directly related to the local transverse correlation 



Tf 1 



D. Properties of weakly coupled ladders 

The interladder coupling J' induces a low temperature 
ordered phase (the 3D-ordered phase in Fig. [2jb). Using 
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FIG. 8. (Color online): Magnetic field dependence of the 
NMR relaxation rate, Tf 1 ^*), at T = 250 mK. The solid 
red line is the bosonization determination using the ladder 
LL parameters for the BPCB couplings (the dashed blue line 
uses the LL parameters of the spin chain mapping) . The black 
circles are the measurements from Ref. I3Q. 



the mean- field approximation presented in Sec. IIIIDl we 
characterize the ordering and compute the critical tem- 
perature and the order parameter related to this phase. 



1. 3D order transition temperature 

In order to compute the critical temperature of the 3D 
transition, we follow Ref. |2| and treat the staggered part 
of the mean-field Hamiltonian Hmf (|2T|) perturbatively 
using linear response. The instability of the resulting 
mean-field susceptibility, due to the 3D transition, ap- 
pears at T c when 71 



xS x (« = o,w = o)| 5 



n c J' 



(26) 



where \a X * s the transverse staggered retarded corre- 
lation function of an isolated single ladder system (ap- 
pendix [B3]). This correlation can be computed analyti- 
cally (see Eq. (|B8|) ) using the LL low energy description 
of the isolated ladder (Eq. ([16]) ) in the gapless regime. 
Applying the condition ([26]) to the LL correlation (|B8|) 
leads to the critical temperature 



FIG. 9. (Color online): Magnetic field dependence of the 
transition temperature between the gapless regime and the 
3D-ordered phase, T c (h z ), is plotted in solid red line for the 
ladder LL parameters of BPCB (in dashed blue line for the 
LL parameters of the spin chain mapping). The NMR mea- 
surements from Ref. [30 are represented by black circles and 
the neutron diffraction measurements from Ref. 31 by green 
dots. 



BPCB (the only free (scaling) parameter in Eq. ([27[) ). 
The asymmetry of the LL parameters induces a strong 
asymmetry of T c with respect to the middle of the 3D 
phase. 

As the mean- field approximation neglects the quantum 
fluctuations between the ladders, the critical temperature 
T c is overestimated for a given J^f- We thus performed 
a Quantum Monte Carlo (QMC) determination of this 
quantity based on the same 3D lattice structure. Let 
us note that QMC simulations of the coupled spin lad- 
der Hamiltonian (pQ) are possible, since the 3D lattice 
structure, Fig. [18j is unfrustrated. In appendix [D] we 
present results on how to determine the critical tempera- 
tures for the 3D ordering transition using QMC. This 
determination shows 31 that the real critical tempera- 
ture is well approximated by the mean-field approxima- 
tion, but with a rescaling of the real interladder coupling 
J' « 27 mK = a -1 J^p witn a ~ °- 74 - Tne rescaling 
factor a is similar to the values obtained for other quasi 
one-dimensional antiferromagnet o 72 i 73 . 



2tt 



f A x J'n c sin {^)B 2 (^,1 
2u 




(27) 



Introducing the computed LL parameters u, K and 
A x (see Fig. [23]) in this expression, we get the criti- 
cal temperature 3 -^ as a function of the magnetic field. 
This is shown in Fig. [9] together with the experimental 
data. This allows us to extract the mean-field interladder 
coupling J^ F w 20 mK for the experimental compound 



2. Zero temperature 3D order parameter 

The staggered order parameter in the 3D-ordered 
phase, m^, can be analytically determined at zero tem- 
perature using the mean-field approximation for the in- 
terladder coupling and the bosonization technique (see 
Sec. lIIID]) . As m x a = y/A^(cos(0(x))) in the bosonization 
description and the expectation valued 3 - of the operator 
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F(K) 



2u 



(28) 



for the sine- Gordon Hamiltonian Hsg ([22]) with 



F(K) 



ir 2 SK 
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(29) 



we can extract 



m x a = ^A X F{K)— 2 I ^ 



(30) 



This can be evaluated in the 3D-ordered phase by intro- 
ducing into ([30]) the LL parameters u, K and A x from 
Fig. [23j Fig. [10] shows the order parameter versus the 
magnetic field determined analytically and numerically 
by DMRG (see Sec. IIIIDlj) . The two curves are al- 
most indistinguishable and exhibit a strongly asymmetric 
camel-like shaped with two maxima close to the critical 
fields. The asymmetry of the curve is again due to the 
presence of the additional triplet states. This asymmetry 
disappears in the spin chain mapping. 



V. DYNAMICAL CORRELATIONS OF A 
SINGLE LADDER 

In the next paragraph, we discuss the possible exci- 
tations and the corresponding correlation functions cre- 
ated by the spin operators. Such dynamical correlations 
allow us to study the excitations of our system. They 
are also directly related to many experimental measure- 
ments (NMR, INS ... ). We first focus on the gapped 
spin liquid and then treat the gapless regime, for fields 
h c i < h z < h C 2 (see Fig. [2]b). All correlations are com- 
puted using the t-DMRG at zero temperature for the 
single ladder (appendix IC 1|) and are compared to ana- 
lytical results when such results exist. In particular we 
check the overlap with the LL description at low energy 
and use a strong coupling expansion (appendix [AJ to 
qualitatively characterize the obtained spectra. We start 
by a discussion of the correlations for the parameters of 
the compound BPCB (see Sec. IVI A|) and then turn to 
the evolution of the spectra with the coupling ratio 7 
from the weak (7 — >> 00) to strong coupling (7 ~ 0). 



A. Zero temperature correlations and excitations 

In a ladder system different types of correlations are 
possible. We focus here on the quantities 



/CO 
dt(SZ qv (t)S§ iqy ) 
-00 



i(ut-ql) 



(31) 




FIG. 10. (Color online): Magnetic field dependence of the 
staggered magnetization per spin, rria(h z ), in the 3D-ordered 
phase determined for the ground state, (a) Its computation 
with the analytical bosonization technique for the LL param- 
eters of the compound BPCB and J' — 27 mK is represented 
by the dash-dotted red line (dashed blue line for the LL pa- 
rameters of the spin chain mapping). The DMRG result for 
J' — 20 mK(= Jmf) is represented by black dots (as a com- 
parison the bosonization result for J' = 20 mK is plotted in 
solid red line). Note, that these two data are almost indis- 
tinguishable, (b) Comparison between NMR measurements 
(black cirles) done at T — 40 mK from Ref. [33 and scaled 
on the theoretical results for J' = 27 mK, neutron diffraction 
measurements on an absolute scale from Ref.[3l| at T — 54 mK 
(T = 75 mK) (red crosses (black dots)) and the computation 
for J' = 27 mK (blue line as in panel (a)). Recent neu- 
tron diffraction measurements as a function of temperature 
suggest that the data of Ref. [3l| was taken at temperatures 
approximately 10 mK higher than the nominal indicated tem- 
perature. 



where Sf = ± S^ 2 are the symmetric (+) and anti- 
symmetric (— ) operators with rung momentum q y = 0, 7r 
and parity P = +1,-1 respectively, a, f3 = 2:,+,—, 
S" q (t) = e lHt Sf 1 e~ lHt (for a single ladder H corre- 
sponds to the Hamiltonian fl2J)) and the momentum q is 
given in reciprocal lattice units a -1 . The rung momen- 
tum q y is a good quantum number. The dynamical cor- 
relations are directly related to INS measurements (see 
Sec. IVI C|) . They select different types of rung excita- 
tions (as summarized in table [I]). Using the reflection and 
translation invariance of an infinite size system (L — >> 00), 
we can rewrite the considered correlations ([3T]) in a more 
explicit form (at zero temperature)^, i.e. 
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TABLE I. Table of the rung excitations created by the sym- 
metric and antisymmetric operators in the decoupled bond 
limit. P is the parity of the operators in the rung direction 
and AM Z the change of the total magnetization. 



where |0) denotes the ground state of H with energy 
Eq, S^ y (q) = e ~ iql ^iq y an d J2\ is the sum over all 
eigenstates | A) of H with energy E\. The form of Eq. ([32|) 
clearly shows that Sq^(q,uj) is non-zero if the operator 
S@ can create an excitation |A) of energy Eq + u and 
momentum q from the ground state. The correlations 
are then direct probes of the excitations |A) in the 
system. 

Since the experimentally relevant case (compound 
BPCB) corresponds to a relatively strong coupling sit- 
uation (7 <C 1, Eq. [8]), we use the decoupled bond limit 
introduced in Sec. Ill CI to present the expected excita- 
tions \t~) or \s). In table [D we summarize the 
rung excitations created by all the operators and their 
properties. For example, the rung parity P is changed by 
applying an operator with rung momentum q y = tt and 
the z-magnetization is modified by AM Z = ±1 by apply- 
ing the operators S^ respectively. 



B. Excitations in the spin liquid 

Using the decoupled bond limit in the spin liquid 
phase, the excitations in the system can be pictured as 
the excitation of rung singlets to rung triplets. At zero 
magnetic field h z = 0, the system is spin rotational sym- 
metric and the different triplet excitations have the same 
energy ~ J±. It has been seen previously that in this 
system both single triplet excitations and two-triplet ex- 
citations play an important rolei^ - — . We discuss these 
excitations in the following focusing on the ones that can 
be created by the symmetric Sq 01 = 2S^ and the an- 
tisymmetric S% a = 2S^ T correlations (see Fig. [11]) for 
the BPCB parameters (|39|) . Note that these correlations 
are independent of the direction a = x,y,z due to the 
spin rotation symmetry. Our results are in very good 
agreement to previous findings^ - —. 
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FIG. 11. (Color online): Momentum-energy dependent cor- 
relation functions at h z = (a) Symmetric part S§ a (q,uj) 
with a — x,y,z. The dashed (black) line marks the (q . uj) 
position of the two-triplet bound state, Eq. (61b) in Ref. fl4l . 
The dash-dotted (white) lines correspond to the boundaries of 
the two-triplet continuum, (b) Antisymmetric part S^ a (q, uj). 
The dashed (black) line corresponds to the predicted disper- 
sion relation of a single triplet excitation (Eq. (8) in Ref. fl3T) . 



1. Single triplet excitations 

At h z = 0, the system is in a global spin singlet stated 
(S = 0). The q y = tt correlation couples then to states 
with an odd number of triplet excitations with rung par- 
ity P = -1 and total spin 5 = 1, M z = ±1,0 (see 
table UJ. Nevertheless, only single triplet excitations are 
numerically resolved. Their spectral weight is concen- 
trated in a very sharp peak whose dispersion relation 
can be approximated using a strong coupling expansion 
in 7. Up to first order it is simply given by a cosine 
dispersion^, i.e. uj t (q)/J± w l±7COsg. Further cor- 
rections up to third order in 7 have been determined in 
Ref.Ojl In Fig.QTJb we compare the numerical results for 
the BPCB parameters (|39j to the expression up to third 
order in 7. The strong coupling expansion describes very 
well the position of the numerically found excitations. 



2. Two-triplet excitations 

The structure of the q y = correlation is more complex 
(Fig. [TTJa). Due to the rung parity P = 1 of the opera- 
tors So, the excitations correspond to an even number of 
triplet excitations with total spin S = 1, M z = ±1,0. 
We focus here on the two-triplet excitations that can 
be resolved numerically. These can be divided into a 
broad continuum and a very sharp triplet (S = 1) bound 
state of a pair of rung triplets. Since these excitations 
stem from the coupling to triplets already present in the 
ground state (Fig. SJ, their amplitude for the considered 
BPCB parameters ([39|) is considerably smaller than the 
weight of the single triplet excitations^. 

The dispersion relation of the bound states has been 
calculated using a linked cluster series expansion 14 . The 



14 



first terms of the expansion have an inverse cosine form 
and the bound state only exists in an interval around 
q = 7r (cf. Ref. [3 and [TBI). The numerical results for 
the BPCB parameters ([39|) agree very well with the an- 
alytic form of the dispersion (Fig. ITTIa). The upper and 
lower limits of the continuum can be determined by con- 
sidering the boundary of the two non-interacting triplet 
continuum. They are numerically computed using the 
single triplet dispersion (Eq. (8) in Ref. fl3l ) and shown 
in Fig. [TTJa. They agree very well with the numerically 
found results. The comparison with the known solutions 
serves as a check of the quality of our numerical results. 

C. Excitations in the gapless regime 

A small applied magnetic field (h z < /i c i), at first or- 
der, only smoothly translates the excitations shown in 
Fig. [TT] by an energy —h z M z due to the Zeeman effect. 
However, if the magnetic field exceeds h c i, the system 
enters into the gapless regime with a continuum of ex- 
citations at low energy. For small values of 7 most fea- 
tures of this low energy continuum are qualitatively well 
described by considering the lowest two modes of the lad- 
der only. Beside the low energy continuum, a complex 
structure of high energy excitations exist. Contrarily to 
the low energy sector this structure crucially depends on 
the high energy triplet modes. In the following, we give 
a simple picture for these excitations starting from the 
decoupled bond limit. 

1. Characterization of the excitations in the decoupled bond 

limit 

The evolution of the spectra for the BPCB parameters 
with increasing magnetic field are presented in Fig.[T2lfor 
S\ z , in Fig. [13] for S+", and in Fig. [H for S"+. Three 
different classes of excitations occur: 

(i) a continuum of excitations at low energy for S§ z 
and 

(ii) single triplet excitations at higher energy with a 
clear substructure for S* z , Sq , and 

(iii) excitations at higher energy for Sq z and Sq~ and 
Sq + stemming from two-triplet excitations which 
have their main weight around q ~ tt. 

In the following we summarize some of the characteristic 
features of these excitations, before we study them in 
more detail in Sees. IVC2l to lVC 3 bl 

(i) The continuum at low energy which does not exist in 
the spin liquid is a characteristic signature of the gapless 
regime. It stems from excitations within the low energy 
band which corresponds to the \s) and states in the 
decoupled bond limit (cf. Fig. Oa and table H]): 

Sq z : excitations within the triplet mode 
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FIG. 12. (Color online): Momentum-energy dependent zz- 
correlation function (1) at m z = 0.25 (h z = 3.153 Jy), 
(2) at m z = 0.5 (h z = 4.194 Jy), and (3) at m z = 0.75 
(h z = 5.192 Jy). (a) Symmetric part So z (q,uj) without Bragg 
peak at q = 0. The dashed black lines correspond to the 
location of the slow divergence at the lower edge of the con- 
tinuum predicted by the LL theory. The dashed white curve 
corresponds to the predicted two-triplet bound state location, 
(b) Antisymmetric part S zz (q, uS). The dashed black lines cor- 
respond to the position of the high energy divergences or cusps 
predicted by the approximate mapping on the t- J model. The 
vertical white dash-dotted lines mark the momenta of the 
minimum energy of the high energy continuum and the black 
cross is the energy of its lower edge 20 at q — tt. 

S^ ± : excitations between the singlet \s) and the triplet 
mode. 

This continuum is smoothly connected to the spin liquid 
spectrum in the case of . It originates from the single 
triplet branch (Sec. IV B 1|) when the latter reaches 
the ground state energy due to the Zeeman effect. Since 
two modes play the main role in the description of these 
low energy features, many of them can already be ex- 
plained qualitatively by the spin chain mapping. The ex- 
citations in the chain have been studied previously using 
a Bethe ansatz description and exact diagonalization cal- 
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FIG. 13. (Color online): Momentum- energy dependent H 

correlation function (1) at m z — 0.25 (h z — 3.153 Jy), 
(2) at m z = 0.5 (h z = 4.194 Jy), and (3) at = 0.75 
(h z — 5.192 J||). (a) Symmetric part Sq (<?,cj). The vertical 
dash-dotted white lines mark the momenta of the minimum 
energy of the high energy continuum and the horizontal ones 
the frequency of its lower edge^ at q = 0, 2tt. The dashed 
white lines correspond to the position of the high energy di- 
vergences or cusps predicted by the approximate mapping on 
the t- J model. The dotted white curve corresponds to the pre- 
dicted two-triplet bound state location. The high energy exci- 
tations at uj = 3h z are hardly visible, (b) Antisymmetric part 
(q,uj). The dashed and dash-dotted (dotted) white lines 
correspond to the location of the strong divergence (cusp) at 
the lower edge of the continuum predicted by the LL theory. 



dilations in Ref. [75|. More recently they were computed 
in Ref. [35| due to recent progress in the Bethe ansatz 
method. In particular, the boundary of the spectrum at 
low energy is well described by this approach, since the 
LL velocity determining it is hardly influenced by the 
higher modes (cf. Fig. [23]). However, a more quantita- 
tive description requires to take into account the higher 
modes of the system as well. In Sec. IV C 21 we compare 
in detail our results with the LL theory and the spin 
chain mapping pointing out their corresponding ranges 



FIG. 14. (Color online): Momentum-energy dependent — h- 
correlation function (1) at m z — 0.25 (h z — 3.153 Jy), 
(2) at m z = 0.5 (h z = 4.194 Jy), and (3) at m z = 0.75 
(h z — 5.192 Jy). (a) Symmetric part (q,uS). The verti- 
cal dash-dotted white lines correspond to the momenta at 
which the minimum energy of the high energy continuum 
occurs and the horizontal line to the frequency of its lower 
edge^ at q — 0, 2tt. The dashed black curve corresponds to 
the predicted two-triplet bound state location, (b) Antisym- 
metric part S~ + (q,uj). The dashed and dash-dotted (dotted) 
white lines correspond to the location of the strong divergence 
(cusp) at the lower edge of the continuum predicted by the 
LL theory. 

of validity. 

(ii) The single high energy triplet excitations form a 
continuum with a clear substructure. In the decoupled 
bond limit, these excitations correspond to 

S zz : Single triplet excitations \t°) at energy ~ h z 

Sq : Single triplet excitations \t°) at energy ~ h z 

S+~ : Single triplet excitations \t~) at energy ~ 2h z . 

Many of the features of these continua can be understood 
by mapping the problem onto a mobile hole in a chain, 
as pointed out first in Ref. j76, We detail in Sec. IV C 3 bl 
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tt(2 - m z ) 



(b) 




7T(1 
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FIG. 15. Fermionic picture for the effect of the magnetic field: 
Filling of (a) the singlet band |s), (b) the triplet band |t + ) in 
the gapless phase for a given magnetization m z . 



and appendix [X] this mapping. It opens the possibility 
to investigate the behavior of a single hole in a t-J like 
model using experiments in pure spin ladder compounds. 

(iii) The high energy continuum, which has almost no 
weight close to the Brillouin zone boundary (q = 0,27r), 
is related to two-triplet excitations of the spin liquid 
(Sec. IV B"2]) . They are generated from high energy triplet 
components of the ground state. Their weight therefore 
vanishes for 7 — )• and the excitations correspond to 

Sq + : Two-triplet excitations. * - \t+)\t )) at 

energy ~ h z 

S§ z : Two-triplet excitations -±={\t + )\t~) - \t~)\t + )) at 
energy ~ 2h z 

: Two-triplet excitations -^(\t°)\t-)-\t-)\t )) at en- 
ergy ~ 3h z . 



2. Low energy continuum 

In this section we concentrate on the low energy exci- 
tations of type (i) discussing first their support and then 
comparing their spectral weight to the LL prediction. 

a. Support of the low energy excitations The po- 
sition of the soft modes in the low energy contin- 
uum can be directly obtained from the bosonization 
representatio n 2 ! 19 i 2Q 



(see appendix IB 2|) . They can also 



be understood in a simple picture which we outline in 
the following. The distribution of the rung state popu- 
lation in the ground state depends on the magnetic field 
h z (see Fig. [4]). Taking a fermionic point of view, the 
magnetic field acts as a chemical potential that fixes the 
occupation of the singlet and triplet rung states. In- 
creasing the magnetic field reduces the number of sin- 
glets, whereas at the same time the number of triplets 
increases (see sketch in Fig. [T5]) . The Fermi level lies 
at the momenta q = 7rm 2; ,7r(2 — m z ) for the singlet 
states and at the momenta q = 7r(l — m z ) 1 ix(l + m z ) 
for the triplet states. In this picture the soft modes cor- 
respond to excitations at the Fermi levels. For transi- 
tions \t + ) <H> \t + ) the transferred momenta of these zero 
energy excitations are q = 0, 2ixm z ^ 27r(l — m z ). In con- 
trast the interspecies transitions \t Jr ) <B> \s) allow the 
transfer of q = 7r(l — 2m 2 ), 7r, 7r(l + 2m z ). Therefore 



the positions of the soft modes in the longitudinal cor- 
relation Sq z which allows transitions within the triplet 
states shift from the boundaries of the Brillouin zone in- 
wards towards q = tt when m z increases (Fig. [12]a). In 
contrast the positions of the soft modes in the transverse 
correlations which allow transitions between the sin- 
glets and the triplets move with increasing magnetic field 
outwards (Figs. [T3lb andfHlb). 

The top of these low energy continua are reached when 
the excitations reach the boundaries of the energy band. 
In particular, the maximum of the higher boundary lies at 
the momentum q = tt which is easily understood within 
the simple picture drawn above (cf. Fig. [T5|) . A more 
detailed description of different parts of these low energy 
continua is given in Ref. [75|. 

Let us compare the above findings with the predictions 
of the LL theory for the dynamical correlation ^ 2 ! 19 ! 20 . 
Details on the LL description of the correlations are 
given in appendix IB 21 The LL theory predicts a lin- 
ear momentum-frequency dependence of the lower con- 
tinuum edges with a slope given by the LL velocity ±u 
(Fig. [23]). The position of the soft modes are given by the 
ones outlined above (see Fig. [23]). The predicted support 
at low energy agrees very well with the numerical results 
(Fig. [12]a, [13]b, and [Hlb). Of course when one reaches 
energies of order Jy in the spectra one cannot rely on 
the LL theory anymore. This is true in particular for the 
upper limit of the spectra. 

b. Spectral weight of the excitations Let us now fo- 
cus on the distribution of the spectral weight in the low 
energy continuum. In particular we compare our numer- 
ical findings to the Luttinger liquid description. Qual- 
itatively, the LL theory predictions for the low energy 
spectra are well reproduced by the DMRG computations. 

The Luttinger liquid predicts typically an algebraic be- 
havior of the correlations at the low energy boundaries 
which can be a divergence or a cusp. 



QZZ 

^0 



: The Luttinger liquid predicts peaks at the q = 0, 2ix 
branches and a slow divergence at the lower edge 
of the incommensurate branches q = 2i\m z ; , 2ix(l — 
m z ) (with exponent 1 — K^0.2< 1). In the nu- 
merical results (Fig. [121 a) a slight increase of the 
weight towards the lower edge of the incommensu- 
rate branches can be seen. 



S+~ : A strong divergence at the lower edge of the q = tt 
branch (with exponent 1 — 1/4K w 3/4 ^> 0) is ob- 
tained within the Luttinger liquid description. This 
is in good qualitative agreement with the strong 
increase of the spectral weight observed in the nu- 
merical data. A more interesting behavior is found 
close to the momenta q = 7r(l db 2m z ) in the in- 
commensurate branches. Here a strong divergence 
is predicted for momenta higher (lower) than the 
soft mode q = 7r(l - 2m z ) (q = 7r(l + 2m z )) 
with exponent 1 — t?_^ 3/4^0. In contrast 
for momenta lower (higher) than the soft mode 
q — 7r(l — 2m z ) (q = 7r(l + 2m z )) a cusp with 
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FIG. 16. (Color online): Cuts at fixed momentum q = tt 
and magnetization m z — 0.5 of the low energy spectrum (a) 
So z (q — ^5^)5 an d (b) ^ + (^ — 7T,u;). The (red) circles and 
the (black) squares are the numerical results for the ladder 
and its spin chain mapping, respectively. The dashed lines 
correspond to the LL predictions and the solid lines are the 
latter convolved with the same Gaussian filter than the nu- 
merical data. The DMRG frequency numerical limitation is 
of the order of the peak broadening of width Sou « 0.1 J\\ (see 
appendix IC 3|) . 

exponent 1 — 77+ ~ — 5/4 <C is expected^. In the 
numerical results (Fig. [13jb) this very different be- 
havior below and above the soft modes is evident. 
The divergence and cusp correspond to a large and 
invisible weight, respectively. 

£~ + : The same behavior as S+~ replacing m z — >> — m z 
can be observed in Fig. [Hlb. 

To compare quantitatively the predictions of the LL 
to the numerical results we show in Fig. [T6l different cuts 
of the correlations at fixed momentum q = tt and mag- 
netization m z = 0.5 for the ladder and the spin chain 
mapping. These plots show the DMRG results, the LL 
description, and the latter convolved with the Gaussian 
filter. The filter had been used in the numerical data 
to avoid effects due to the finite time-interval simulated 
(see appendix IC 3|) . Note that the amplitude of the LL 
results are inferred from the static correlation functions, 
such that the LL curve is fully determined and no fitting 
parameter is left. Therefore the convolved LL results 
can directly be compared to the numerical results. Even 
though the presented numerical resolution might not be 
good enough to resolve the behavior close to the diver- 
gences (cusps), interesting information as the arising dif- 
ferences between the spin chain mapping and the full lad- 



der calculations can already be extracted. In Fig. [T6la, 
we show a cut through the correlation at fixed momen- 
tum S zz (q = 7r,co>). The convolved LL and the numerical 
results compare very well. The difference between the 
real ladder calculations and the spin chain mapping that 
neglect the effects of the higher triplet states \t°) is 
obvious. From the LL description point of view, the shift 
of the spin chain correlation compared to the real ladder 
curve comes mainly from the prefactor A z and the alge- 
braic exponent which are clearly modified by the effects 
of the high energy triplets (see Fig. |23|) . 

For the transverse correlations, the LL theory predicts 
a strong divergence (with an exponent 1 — 1 / AK w 3/4 ^> 
at the lower boundary of the continuum branch at q = 
7r. A cut through the low energy continuum S~ + (q = 
7t,cl>) is shown in Fig. [16jb. The convolved Luttinger 
liquid reproduces well the numerical results. 



3. High energy excitations 

a. Weak coupling description of the high energy ex- 
citations Before looking in detail at the two kinds of 
high energy excitations presented in Sec. IV C II we com- 
pare our computed high energy spectra with the weak 
coupling description (7 ^> 1). In this limit information 
on the spectrum can be extracted from the bosonization 
description ^ 19 i 2Q . In particular one expects a power law 
singularity at the lower edge continuum with a minimal 
position at q = 7r(l ± m z ) (for S zz ), q = 7rm z , 7r(2 — m z ) 
(for S^ T ) and an energy h z at momentum q = tt (for 
S zz ), q = (for S^ T ). Except for Sq + in which the spec- 
tral weight is too low for a good visualization, our com- 
puted spectra reproduces well the predictions for the min- 
imal positions even though the coupling strength consid- 
ered is not in the weak coupling limit (cf. Figs.[T2lb.[T3la 
and El a). 

b. High energy single triplet excitations The high 
energy single triplet continua originate from the tran- 
sition of the low energy rung states \s) and to the 
high energy triplets \t°) and \t~). The excitations com- 
ing from the singlets \s) (in S zz and S+~) are already 
present in the spin liquid phase (cf. Sec. IV B 1)) in which 
they have the shape of a sharp peak centered on the 
triplet dispersion. The transition between the gapped 
spin liquid and the gapless regime is smooth and consists 
in a splitting and a broadening of the triplet branch that 
generates a broad continuum of new excitations. Con- 
trarily to the latter the excitations coming from the low 
energy triplets (in Sq ) are not present in the spin 
liquid phase. The corresponding spectral weight appears 
when h z > h c \. 

An interpretation of the complex structure of these 
high energy continua can be obtained in terms of itiner- 
ant quantum chains. Using a strong coupling expansion 
of the Hamiltonian (|2]) (appendix [A]) it is possible to map 
the high energy single triplet excitations \t°) to a single 
hole in a system populated by two types of particles with 
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pseudo spin |f) = ||) = |s) (with the Sec. Ill CI no- 
tation). 

In this picture the effective Hamiltonian of the J± 
energy sector is approximately equivalent to the half 
filled anisotropic ID t-J model with one hole (see ap- 
pendix [A3T]) . The effective Hamiltonian is given by 



Hi 



t-j 



H 



xxz 



(33) 



where e = (Ji + /i 2 )/2 is an energy shift and H t = 

J\\/2 J2i a ( c l a c i+i a + ^i.c.) is the usual hopping term. 

Here cj a (q j<t ) is the creation (annihilation) operator of 

a fermion with pseudo spin a = f,| at the site I. Note 
that although we are dealing here with spin states, it is 
possible to faithfully represent the three states of each 
site's Hilbert space (|s), \t )) using a fermion rep- 

resentation. 

Additionally to the usual terms of the t-J model a near- 
est neighbor interaction term between one of the spins 
and the hole arises 



H s - h = ~ ^2 [ n ^ n z+i,t + n iA n i+^h 



(34) 



Here n^h is the density operator of the hole on the site 
I. In this language the spectral weight of S zz and Sq~ 
corresponding to the single high energy triplet excitations 
is equivalent to the single particle spectral functions of 
the up-spin and down-spin particle respectively: 

S zz oc (ctc~) with hole of type \s) \t°) 

S$~ oc (c|c-> with hole of type \t+)^\t°). ^ 

Here <40(^) = Ea \(XKa\0)\ 2 S(uj + E - E x ). 

For the standard t-J model (for SU (2) invariant XXX 
spin background and without the anisotropic term H s _h 
in Eq. (j33j), these spectral functions have been studied 
in Refs. M and The presence of singularities of the 
form 



(40fe^) oc [uj-uj t o(q - q u )] 



2X v (q)-l 



(36) 



were found. Here uj t o(q) is the \t°) triplet dispersion re- 
lation, q v the spinon momentum at the Fermi level and 
X v the algebraic decay exponent at the singularity. This 
exponent is not known in our case and depends on the 
magnetization m z and the momentum q. It generates 
a peak or a cusp at the energy u = u t o(q — q v ). The 
spinon momentum q v depends on the type of the rung 
state before excitation (y = s, t + ). For an excitation cre- 
ated from a singlet state q s = ±7rm z (for S zz ) and from 
the triplet state q t + = tt(1 ± m z ) (for (Fig. [15]) . 

At h z = 0, a series expansion of u t o(q) can be performed 
(Eq. 8 in Ref. 13). To extend it into the gapless phase 
(h c i < h z < h C 2), we approximate cj £ o(g) by shifting the 
value LJt(q) at h z = by the Zeeman shift, i.e. 



uj t o(q) =oj t (q) + AE (h z ). 



(37) 



Here we used the shift of the ground state energy per 
rung AE (h z ) = E (h z ) - E (0). AE was determined 
by DMRG calculations (Fig. [251b for the BPCB pa- 
rameters). The resulting momentum- frequency positions 
uj = uj t o(q — qv) of the high energy singularities (cusps or 
divergencies) are plotted on the spectrum Fig. [12jb and 
Fig. [13] a. They agree remarkably well with the shape 
of the computed spectra in particular for small magnetic 
field 80 . Neglecting the additional interaction term 
the t-J model Hamiltonian would lead to a symmetry 
of these excitations with respect to half magnetization. 
However in the numerical spectra the effect of the inter- 
action shows up in a clear asymmetry of these excitations 
(Figs. [121 l.b and [T3l3.a). In particular in the 5q cor- 
relation some of the weight is seemingly detaching and 
pushed towards the upper boundary of the continuum 
(Fig. [13l3.a) for large magnetization. A more detailed 
account of the spectra can be found in Ref. [8l|. 

A similar mapping can be performed for the single \t~) 
excitation. In contrast to the J± sector in the 2J± sec- 
tor not only the excitation occurs, but the effective 
Hamiltonian mixes also \t°) triplets into the description. 
Therefore the description by a single hole in a spin- 1/2 
chain breaks down and more local degrees of freedom 
are required. This results in a more complex structure 
as seen in Fig. [13jl.b. Previously high-energy excita- 
tions in dimerized antiferromagnets have been described 
rather generally by a mapping to an X-ray edge singu- 
larity problem^"—. It is interesting though that in the 
present setup these excitations can be understood as t- 
J hole spectral functions, which display a much richer 
structure than anticipated. 

c. High energy two-triplet excitations The two- 
triplet continua and bound states already discussed in 
the spin liquid phase (cf. Sec. IV B 2|) are still visible in 
the gapless regime in the symmetric correlations (S zz and 
S^ T ). At low magnetic field the location of their maxi- 
mal spectral weight can be approximated by the expres- 
sion of the bound state dispersion at zero field (Eq. (61b) 
in Ref. Il4[ ) shifted by the Zeeman energy 85 . The two- 
triplet excitation location obtained in this way agrees to a 
good extent with the location found in the numerical cal- 
culations (cf. Figs.[T2ll.a.[T3ll.a and[T4ll.a). Since these 
excitations are generated from the high energy triplet 
components in the ground state and these vanish with in- 
creasing magnetic field (cf. Fig. 2]), their residual spectral 
weight slowly disappears with increasing magnetization. 



D. Weak to strong coupling evolution 

For all the excitation spectra presented above the intra- 
chain coupling ratio of BPCB 7 = J\\/J± ~ 1/3.55 <C 1 
was taken. For this chosen value of 7, a strong coupling 
approach gives a reasonable description of the physics. In 
this section we discuss the evolution of the spectra from 
weak (7 00) to strong coupling (7 — >• 0). To illustrate 
this behavior, we show in Fig. [TTlthe symmetric and an- 
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tisymmetric parts of the correlations 5+ at m z = 0.25 
for different coupling ratios 7 = oo,2,l,0.5,0. 

At 7 — > oo (Fig. [13 1), the chains forming the ladder 
correspond to two decoupled Heisenberg chains. In this 
case the symmetric and antisymmetric correlations are 
identical Sq~ = S+~ and are equivalent to the correla- 
tion 25^ of the single chain^ with magnetization per 
spin m z /2 = 0.125. A complex low energy continuum 
exists with zero energy branches ^ 19 i 75 at momenta q = 
±7rm z ,7r similar to that discussed in Sec. IV C 21 In con- 
trast, in the strong coupling limit (7 — >> 0) (Fig. [T7l5.b) 
the symmetric correlations vanish and the antisymmet- 
ric part corresponds to the single chain correlation 2S^ 
with anisotropy A = 1/2 and magnetization per spin 
m z — 1/2 (see the spin chain mapping in Sec. Ill C|) . The 
antisymmetric part consists of a low energy continuum 
with branches at momenta q = (l±2m 2 )7r, tt (Sec. IV C 2|) . 
Note, that a bosonization description of the low energy 
sectors of both extreme regimes can be formulated ^ 19 i 2Q 
(appendix IB 2ft . 

In the following we discuss the evolution between 
these two limits. In the antisymmetric correlation 
(cf. Fig. [171 2-5. b) a low energy continuum exists at all 
couplings with a zero energy excitation branch at q = tt. 
These low energy excitations correspond mainly to the 
excitations with AS = AM Z = —1. This has been 
pointed out for the weak coupling limi t 19 ! 75 . They be- 
come the transitions — >• \s) with the same quan- 
tum numbers in the decoupled bond limit. Addition- 
ally the upper part of the excitation spectrum at weak 
coupling, which mainly corresponds to excitations with 75 
AS = 0, 1 and AM Z = — 1 splits from the lower part of 
the spectrum and moves to higher energy while increas- 
ing the coupling. It evolves to a high energy excitation 
branch which corresponds in the decoupled bond limit to 
the 1 5) — >• \t~) transition, i.e. single triplet excitations of 
type (ii) fSec. lVC3b|) approximately at 86 2J ± . 

The properties of the zero energy excitation branch 
at q = tt show a smooth transition between the two 
limits 2,19 . For example the slope of the lower edge con- 
tinuum which is determined by the LL velocity u de- 
creases smoothly from its value for the Heisenberg chain 
to the lower value for the anisotropic spin chain with 
A = 1/2 in the strong coupling limit. In contrast to 
this smooth change, the presence of a finite value of J± 
leads to the formation of a gap in the incommensurate 
low energy branches^ at q = ±.irm z . With increasing 
coupling strength J± new low energy branches at mo- 
menta q = 7r(l ± 2m z ) become visible 2,20 . The weight 
of these gapless branches is very small for small coupling 
and increases with stronger coupling 2 -. 

In contrast to the antisymmetric part, the symmet- 
ric part Sq becomes gapped when the interladder cou- 
pling J± is turned on. The lowest energy excitations 
remain close to the momenta q = ±.irm z in agreement 
with Ref. l20|. They connect to the single triplet excita- 
tions of type (ii) (Sec. IV C 3b[) which are approximately 
at an energy J±. While increasing 7 the higher part of 



the spectrum starts to separate from the main part and 
evolves to a branch of high energy two-triplet excitations 
of type (hi) (Sec. I V C 3 c|) . These are located at approx- 
imately 3J±. Our computed spectra for 7 = 2,1,0.5 
presented in Fig. [1712-4. a clearly show this behavior. In 
Fig. [T7l4.a the highest two-triplet excitations cannot be 
seen anymore since their spectral weight is too low. 



E. Influence of the weak interladder coupling on 
the excitation spectrum 

Up to now we only discussed the excitations of a single 
spin ladder and neglected the weak interladder coupling 
J' usually present in real compounds. 

Deep inside of the spin liquid phase, the correlations 
for a single ladder are dominated by high energy single 
or multi triplet excitations as discussed in Sec. IV Bl The 
presence of a small interladder coupling J' causes a dis- 
persion in the interladder direction with an amplitude of 
order J 7 . This effect can be evaluated for independent 
triplet excitations using a single mode approximation^. 
However, for the compound BPCB the interladder cou- 
pling is so small that present day experiments do not 
resolve this small broadenin g 32 ! 33 . 

In contrast in the gapless phase the effect of the in- 
terladder coupling can change considerably the excita- 
tions. In particular below the transition temperature to 
the 3D-ordered phase, a Bragg peak appears at q = tt in 
the transverse dynamical functions S^ T . As discussed in 
Ref. [59), this Bragg peak is surrounded by gapless Gold- 
stone modes and it has been measured in the compound 
BPCB 31 . Additional high energy modes are predicted 
to occur in the transverse S^ T and longitudinal Sq z59 . 
It would be interesting to compute the excitations using 
random phase approximation analogously to Ref. [59| in 
combination with the computed dynamical correlations 
for the single ladder in order to investigate the effect of a 
weak interladder coupling in more detail. However, this 
goes beyond the scope of the present work and will be 
left for a future study. 



VI. EXPERIMENTAL MEASUREMENTS 

In this section we summarize the results of experiments 
on the compound BPCB and compare them to the the- 
oretical predictions. First we introduce the structure of 
the BCPB compound. Then we focus on static and low 
energy results of the system. Finally we discuss the de- 
tection of excitations by INS 32 i 33 measurements. In par- 
ticular we give a prediction of the INS cross section in 
a full range of energy and compare its low energy part 
with the measured spectra on the compound BPCB. 
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FIG. 17. (Color online): Momentum- energy dependent H correlations (S^ y ~(q, lj)) at m z = 0.25 for different ladder couplings 

(7 = J\\/J±) (1) 7 — >• 00, (2) 7 = 2, (3) 7 = 1, (4) 7 = 0.5, (5) 7 — >• 0. The symmetric (antisymmetric) correlations with q y = 
(q y = 7r) are presented in the figures labeled by a (b). In (1,2-4. a) the vertical dashed lines represent the incommensurate 
momenta of the low energy branches of the single spin chain at q = ±Tim z — ±7r/4 (they also correspond to the predicted 
momenta of the lowest energy excitations of the symmetric correlations^). The horizontal solid (dotted) horizontal lines in 
(2-4. a) correspond to the approximate energy J± (3Jj_) of the single triplet excitations of type (ii) (two-triplet excitations of 
type (hi)). The horizontal dashed lines in (2-4. b) correspond to the approximate energy 2J± of the excitations of type (ii). The 
dash-dotted lines in (1) correspond to the linear low energy boundaries of the continuum of excitations given by the LL theory 
applied on single Heisenberg chains 7 — >> 00. The dash-dotted lines in (2-5. b) correspond to the linear low energy boundaries 
of the continuum of excitations given by the LL theory on spin ladder with finite 7 ( appendix IB 2|) . 



A. Structure of the compound (C 5 Hi 2 N) 2 CuBr4 

The compound (C 5 Hi2N) 2 CuBr4, customarily 
called BPCB or (Hpip)2CuBr4, has been inten- 
sively investigated using different experimental 
methods such as nuclear magnetic resonance 3 -^ 
(NMR), neutron diffraction^ (ND), inelastic neutron 
scatterin g 32 ! 33 (INS), calorimetry 2 ^, magnetometry^i, 
magnetostrictio n 88 ! 89 , and electron spin resonance 
spectroscopy 9 ^ (ESR). 

The magnetic properties of the compound are related 
to the unpaired highest energy orbital of the Cu 2+ ions. 
Thus the corresponding spin structure (Fig. [T8j) matches 
with the Cu 2+ locatio n 28 ! 30 ! 32 . The unpaired spins form 
two types of inequivalent ladders (Fig. [T8|) along the a 
axis (a, b and c are the unit cell vectors of BPCB). 
The direction of the rung vectors of these ladders are 
di, 2 = (0.3904, ±0.1598, 0.4842) in the primitive vector 
coordinates (Fig. fl8lb). As one can see from the projec- 
tion of the spin structure onto the bc-plane (Fig. H51b), 
each rung has n c = 4 inter ladder neighboring spins. 

The BPCB structure has been identified as a good ex- 
perimental realization of the system of weakly coupled 
spin- 1/2 ladders 2 ^ described by the Hamiltonian (pQ). As 



we will explain in the next section, the interladder cou- 
pling J' has been experimentally determined to b o 3Q i 31 

J' « 20 - 100 mK. (38) 

The intraladder couplings from Eq. (J2j) were determined 
to be J ± w 12.6 - 13.3 K, J y w 3.3 - 3.8 K with dif- 
ferent experimental techniques and at different experi- 
mental conditions 2 ^ -33 1 87 ~ — . In this paper, we use the 
values 91 

J± « 12.6 K, J\\ « 3.55 K. (39) 

Recently a slight anisotropy of the order of 5% of 
J_l has been discovered by ESR 9 -^ measurements. This 
anisotropy could explain the small discrepancies between 
the couplings found in different experiments. The mag- 
netic field in Tesla is related to h z replacing 

h z -+ g^ B h z (40) 

in Eq. ([2]) with \±b being the Bohr magneton and g being 
the Lande factor of the unpaired copper electron spins. 
The latter depends on the orientation of the sample with 
respect to the magnetic field. For the orientation chosen 
in the NMR measurements 3 ^, it amounts to2i g « 2.126. 
It can vary up to ~ 10% for other experimental setups 2 ^. 
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FIG. 18. (Color online): Coupling structure of BPCB where 
the unpaired electron spins of the Cu 2+ atoms in the first 
(second) type of ladders are pictured by red (blue) spheres. 
The J_l, J\\ and J' coupling paths are represented in turquoise, 
pink, and green, respectively, a, b and c are the three unit cell 
vectors of the structure. Gray arrows are the rung vectors of 
the two types of ladders di,2. (a) 3D structure, (b) Projection 
of the 3D structure onto the bc-plane. (c) Projection of the 
3D structure onto the ac-plane. (d) Projection of the 3D 
structure onto the ab-plane. 



B. Thermodynamic measurements and low energy 
properties 

Many interesting thermodynamic measurements have 
been performed on BPCB. We select in the following 
some of these experiments and compare them to the the- 
oretical predictions. 

The longitudinal magnetization that can be measured 
very precisely by NMR (see Ref. |30| ) was shown to agree 
remarkably with the one computed using the weakly cou- 
pled ladder model (Fig. [3)). Unfortunately, the magne- 
tization is not very sensitive to the underlying model 
(Sec. IIV A|) . Thus it cannot be used to select between 
various models, but once the model is fixed, it can be 
used to fix precisely the parameters given the high accu- 
racy of the experimental data. In particular, the position 
of the critical fields are very sensitive to the values of the 
intraladder couplings (Sec. IIV X)) . The couplings deter- 
mined by this method are J± « 12.6 K and Jy ~ 3.55 K. 

A more selective test to distinguish between various 
models is provided by the specific heat. This is due to 
the fact that the specific heat contains information on 
high energy excitations which are characteristic for the 
underlying model. As shown in Fig. [6] the experimen- 
tal data are remarkably described, up to an accuracy of 
a few percent, by a simple Heisenberg ladder Hamilto- 
nian with the parameters extracted from the magneti- 



zation. In particular, not only the low temperature be- 
havior is covered by the ladder description, but also the 
higher maxima. This indicates that the ladder Hamilto- 
nian is an adequate description of the compound. The 
small discrepancies between the specific heat data and 
the calculation which is essentially exact can have var- 
ious sources. First of all, the substraction of the non- 
magnetic term in the experimental data can account for 
some of the deviations. Furthermore the inter ladder cou- 
pling can induce slight changes in the behavior of the 
specific heat. Finally deviations from the simple lad- 
der Hamiltonian can be present. Small anisotropy of the 
couplings can exist and indeed are necessary to interpret 
recent ESR experiments 90 . Other terms such as longer 
range exchanges or Dzyaloshinskii-Moryia (DM) terms 
might occur along the legs even if the latter is forbidden 
by symmetry along the dominant rung coupling. Clearly 
all these deviations from the Heisenberg model cannot be 
larger than a few per cents. They will not lead to any size- 
able deviation for the Luttinger parameters (Fig. [23]) in 
the one dimensional regime. Close to the critical points 
they can, however, play a more important role. It would 
thus be interesting in subsequent studies to refine the 
model to take such deviations into account. 

After having fixed the model and the intraladder cou- 
plings up to a few percents we use it to compute other 
experimentally accessible quantities such as the magne- 
tostriction, thermal expansion, the NMR relaxation rate, 
the transition temperature to the ordered phase and its 
order parameter. In Refs . l88l and l89l the magnetostriction 
and thermal expansion were compared to the theoretical 
results using the described ladder model. A very good 
agreement was found in a broad range of temperature 
(not shown here). Note that only a full ladder model 
allows a global quantitative description of the magne- 
tostriction effect which provides an additional confirma- 
tion of the applicability of the model. The quality of the 
determination of the model and its intraladder param- 
eters becomes even more evident in the comparison of 
the NMR data for the relaxation rate T x _1 with the the- 
oretical results of the Luttinger liquid theory as shown 
in Fig. [HJ Only one adjustable parameter is left, namely 
the hyperfine coupling constant (see Sec. IIV C]) . This pa- 
rameter allows for a global expansion of the theoretical 
curve, but not for a change of its shape. The agreement 
between the theory and the experimental data is very 
good over the whole range of the magnetic field and only 
small deviations can be seen. The compound BPCB thus 
allows to quantitatively test the Luttinger liquid univer- 
sality class. Even though the Luttinger liquid description 
is restricted to low energies, in BPCB its range of validity 
is rather large. Indeed at high energy, its breakdown is 
approximately signaled by the first peak of the specific 
heat 29 (see Sec. I IVB 2[) . Here this has a maximum scale 
of about T ~ 1.5 K at midpoint between h& and h C 2 (see 
Fig. 0c). Given the low ordering temperature which has 
a maximum at about T ~ 100 mK this leaves a rather 
large Luttinger regime for this compound. 
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Taking now the coupling between ladders into account, 
one can induce a transition to a three-dimensional or- 
dered phase. The transition temperature is shown in 
Fig. [9l Experimentally it is determined by NMR££ and 
neutron diffraction measurements 3 -^. Theoretically the 
ladders are described by Luttinger liquid theory and their 
interladder coupling is treated in a mean- field approxima- 
tion fSecs. lIIICl and llllDj) . As shown in Fig. El the Lut- 
tinger liquid theory provides a remarkable description of 
the transition to the transverse antiferromagnetic order 
at low temperatures. The shape of T c (h z ) is almost per- 
fectly reproduced, in agreement with both the NMR 30 
and the ND data 3 ^. The comparison with the exper- 
iments determines the interladder coupling J', the only 
adjustable parameter. The simple mean- field approxima- 
tion would give a value of J' ~ 20 mK. As discussed in 
Sec. IIVD II mean-field tends to underestimate the cou- 
pling and it should be corrected by an essentially field 
independent factor. Taking this into account we obtain 
a coupling of the order of J' = 27 mK. 

The order parameter in the antiferromagnetic phase 
can also be observed by experiments. It shows a very 
interesting shape. At a pure experimental level neutron 
diffraction and NMR have some discrepancies as shown 
in Fig. [TUJ These discrepancies can be attributed to the 
different temperatures at which the data has been taken, 
and a probable underestimate of the temperature in the 
neutron diffraction experiments 3 -^. Indeed the order pa- 
rameter close to the critical magnetic field h C 2 is very 
sensitive to temperature, since the transition tempera- 
ture drops steeply in this regime. Note that although 
the NMR allows clearly for a more precise measurement 
of the transverse staggered magnetization it cannot give 
its absolute value. Thus the amplitude of the order pa- 
rameter is fixed from the neutron diffraction measure- 
ment. Even though a good agreement between the the- 
oretical results and the experimental results is obtained, 
several questions concerning the deviations remain to be 
addressed. 

First, the theoretical curve does not fully follow the 
shape of the experimental data. Particularly at high 
fields the experimental data shows a stronger decrease. 
A simple explanation for this effect most likely comes 
from the fact that the calculation is performed at zero 
temperature, while the measurement is done at 40 mK. 
This is not a negligible temperature with respect to T c , in 
particular at magnetic fields close to h C 2. Extrapolation 
of the experimental data to zero temperature^ improves 
the agreement. However, for a detailed comparison either 
lower temperature measurements or a calculation of the 
transverse staggered magnetization at finite temperature 
would be required. Both are quite difficult to perform 
and will clearly require further studies. 

The second question comes from the amplitude of the 
staggered magnetization. Indeed the experimental data 
seem to be slightly above the theoretical curve, even if one 
uses the value J' = 27 mK for the interladder coupling. 
This is surprising since one would expect that going be- 



yond the mean- field approximation could only reduce the 
order parameter. Naively, one would thus need a larger 
coupling, perhaps of the order of J' ~ 60 — 80 mK to 
explain the amplitude of the order parameter. This is a 
much larger value than the one extracted from the com- 
parison of T c . How to reconcile these two values remains 
open. The various anisotropics and additional small per- 
turbations in the ladder Hamiltonian could resolve part 
of this discrepancy. However, it seems unlikely that they 
result in a correction of J' by a factor of about 2-3. An- 
other origin might be the presence of some level of frus- 
tration present in the interladder coupling. Clearly more 
experimental and theoretical studies are needed on that 
point. 



C. Inelastic neutron scattering 

The inelastic neutron scattering (INS) technique is a 
direct probe for dynamical spin-spin correlation func- 
tions. Measurements have been performed on the com- 
pound BPCB in the spin liquid phas o 32 i 33 (low magnetic 
field) and in the gapless regime 3 - 2 -. Modeling the com- 
pound BPCB by two inequivalent uncoupled ladders ori- 
ented along the two rung vectors di^ (see Fig. [T8|) the 
magnetic INS cross section^ 2 - is given by the formula 



d 2 a 
dttdE 



[cos(Q • di)(S+- + 5 2 -+) + S+f 
[cos(Q • dOSff 



+2 1 



■5ff]}. (41) 



Here Q = (Q x , Q y ,Q z ) = q — q' is the momentum trans- 
ferred to the sample (q, q' are the incident, scattered 
neutron momenta) and uo = E — E r is the transferred en- 
ergy (E, E' are the incident, scattered neutron energies). 
The correlations Sff are understood to be evaluated at 
the momentum Q • a and frequency cj, and are defined 

by 



S?f (Q-a )W ) = 



Sf(Q-a,w)±S^(Q-a,cu) 



(42) 



(+ sign if i = j and — sign if i ^ j) with defined at 
zero temperature in Eq. ([3T]) and evaluated for a momen- 
tum q = Q • a along the a unit cell vector (momentum 
along the ladder direction). The magnetic form factor 
F(Q) of the Cu 2+ and the ratio q'/q are corrected in the 
experimental data. 

The INS cross section (|4Tj) is directly related to a 
combination of different correlation functions S^^ with 
weights depending on the transferred momentum Q and 
the magnetic field orientation. In the model definition 
(see Sec. Ill Ap . the magnetic field h is pointing along the 
z direction. Aligning it to the b unit cell vector and 
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FIG. 19. (Color online): Theoretical momentum-energy de- 
pendent INS cross section for BPCB with Q = tt (i — 1,2) 
and q — Q • a at m z = 0.5 in (a) a ladder system and (b) the 
spin chain mapping. The horizontal dashed lines correspond 
to the constant energy scans at uj — 0.2,0.4 meV shown in 

Fig. eh 

tuning Q in the a*c*-plane (a*, b* and c* are the recip- 
rocal vectors of a, b and c) allows one to keep constant 
the prefactors in Eq. (|4T|) scanning the a*-momenta with 
the condition Q • di = or tt to target the symmetric or 
antisymmetric part, respectively. 

We focus here on the antisymmetric part for which 
the low energy spectra have already been studied exper- 
imentally and theoretically 32 . Theoretically the focus so 
far lay on the description by the spin chain mapping. We 
compute here the INS cross section (|4T|) for the full ladder 
at m z = 0.25, 0.5, 0.75 using the correlations presented in 
Sec. El The results are shown in Figs. fT9l and [20l and are 
compared to the results from the spin chain mapping. 

As expected from expression (|4T]h it contains the dif- 
ferent excitations present in the spectra of S zz and 
(cf. Figs. H2b, lEflb and [Hlb): 

(a) The low energy continuum originates from the 
transversal correlations . It is qualitatively 
well described by the spin chain mapping that 
presents a symmetry with respect to half magne- 
tization. 

(b) The continuum of excitations at energy 86 ~ J±_ 
comes from the longitudinal correlation S zz and is 
not present in the spin chain mapping. 

(c) The continuum of excitations at energy 86 ~ 2J± 
comes from the transversal correlation and is 
not present in the spin chain mapping. 

The main features of the low energy continuum (a) 
are well covered by the spin chain mapping 32 . However, 
slight differences between the low energy excitations in 
the spin ladder and the spin chain are still visible (cf. also 
Sec. IV C 2|) . They can even be distinguished in the ex- 
perimental data as shown in Figs. [2T1 and [22l where some 
cuts at fixed energy uj = 0.2,0.4meV are plotted. The 



INS measured intensity is directly compared to the the- 
oretical cross section computed for the ladder and 
the spin chain mapping at m z = 0.24, 0.5, 0.72 convolved 
with the instrumental resolution. The amplitude is fixed 
by fitting one proportionality constant for all fields, en- 
ergies, and wave vectors. 

These scans at fixed energy present peaks when 
the lower edge of the continua (related to the cor- 
relations S^ T ) is crossed (see dashed white lines in 
Figs. IT9l and [20]). As one can see, the theoretical curves 
for the ladder and the spin chain both reproduce well the 
main features in the experimental data and only small 
differences are present: 

- The spectral weight intensity at m z = 0.5 and uj = 
0.4 meV (in Fig. [2TJb) is slightly overestimated by 
the spin chain mapping. 

- The height of the two central peaks at m z =0.24 
and uj = 0.2 meV (in Fig. [22jc) is underestimated 
by the spin chain mapping. 

Whereas the low energy excitations (a) only showed a 
slight asymmetry with respect to the magnetization, a 
very different behavior can be seen in the high energy 
part (b)-(c). Indeed, the high energy part of the INS 
cross section (Fig. [20]) is very asymmetric with respect 
to half magnetization. As we discussed in Sec. El these 
excitations are due to the high energy triplets that can be 
excited in S zz and S+~ (see Fig. [12jb and Fig. [T51b) and 
are totally neglected in the spin chain mapping. Their 
corresponding spectral weight is of the same order than 
the low energy spectra, and thus should be accessible 
experimentally. It would be very interesting to have an 
experimental determination of this part of the spectrum, 
since as we have seen it contains characteristic informa- 
tion on the system itself and related to itinerant systems 
via the various mappings. 



VII. CONCLUSIONS 

In this paper we have looked at the thermodynamic 
and dynamical properties of weakly coupled spin lad- 
ders under a magnetic field. This was done by a com- 
bination of analytical techniques, such as Bet he ansatz, 
bosonization and Luttinger liquid theory, and numerical 
techniques such as DMRG and quantum Monte Carlo. 
Using this combination of techniques we were able to 
explore the physical properties in the three main re- 
gions of the phase diagram of such spin ladders under 
a magnetic field, namely: a) a gapped spin liquid at 
low fields; b) a massless phase at intermediate fields 
h c \ < h z < h C 2\ c) a saturated phase at larger fields. 
In addition to the theoretical analysis we compared our 
findings to the experimental results on the compound 
BPCB ((C 5 Hi2N) 2 CuBr4) which is an excellent realiza- 
tion of such ladder systems. 
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FIG. 20. (Color online): Theoretical momentum- energy dependent INS cross section for BPCB with Q • di = ir (i = 1, 2) and 
q = Q a, (a) at m z = 0.25 (b) at m z = 0.25, 0.75 in the spin chain mapping, (c) at m z = 0.75. The horizontal dashed lines 
correspond to the constant energy scans at uj — 0.2, 0.4 meV plotted in Fig. 1221 




FIG. 21. (Color online): Inelastic neutron scattering intensity 
measured along a* of BPCB^ with the momentum tt in the 
rung direction (Q • di = tt) at h z = 10.1 T (m z « 0.5) and 
T — 250 mK after substraction of the zero- field background. 
In each panel, fixed energy scans (shown by white dashed lines 
in Fig. HH) are plotted: (a) u = 0.2 meV, (b) u = 0.4 meV. 
The circles correspond to the experimental data. The red 
(black) solid lines are the m z = 0.5 theoretical data for the 
ladder (the spin chain mapping) convolved with the instru- 
mental resolution. The shaded bands indicate the error bar 
in the experimental determination of a single proportionality 
constant valid for all fields, energies, and wave vectors. The 
width of these areas combines the statistics of all our scans 
with uncertainties in the exact magnetization values at the 
chosen fields and in the convolution procedure. 



For thermodynamics we computed the magnetization 
and specific heat of the system as a function of tempera- 
ture and magnetic field. The extension of the DMRG 
technique to finite temperature allows us to compute 
these quantities with an excellent accuracy. In the gap- 
less phase the low energy part of the specific heat agrees 
well with the prediction of the Luttinger liquid theory, 
which is the low energy theory describing most of mass- 
less one dimensional systems. At higher temperatures, 
the numerical solution is needed to capture the precise 
structure of the peaks in the specific heat, that reflect 



the presence of the excited states in the ladder. Com- 
parison of the theoretical calculations with the measured 
magnetization and specific heat proves to be remarkable. 
This good agreement confirms that the ladder model is 
indeed a faithful description of the compound BPCB. It 
also gives direct access, via the maxima in magnetization 
and peaks in the specific heat to the approximate region 
of applicability of the Luttinger liquid description. For 
the low energy dynamics we used a combination of the 
numerical techniques to determine the Luttinger liquid 
parameters and then the analytical description based on 
Luttinger liquids to compute the dynamical spin-spin cor- 
relation functions. This allowed to extract the NMR re- 
laxation rate and the one dimensional antiferromagnetic 
transverse susceptibility. If the ladders are very weakly 
coupled, which is the case in the considered material, the 
divergence of the susceptibility leads to a three dimen- 
sional antiferromagnetic order at low temperatures in the 
direction transverse to the applied magnetic field. We 
computed this transition temperature and the order pa- 
rameter at zero temperature. Comparison with the mea- 
sured experimental quantities both by NMR and neutron 
diffraction proved again to be remarkable. This excel- 
lent agreement between theory and experiment for these 
quantities as a function of the magnetic field allows to 
quantitatively test the Luttinger liquid theory. Indeed it 
shows that several different correlations are indeed totally 
described by the knowledge of the two Luttinger liquid 
parameters (and the amplitudes relating the microscopic 
operators to the field theory one). This is something that 
could not really be tested previously since either the mi- 
croscopic interactions were not known in details, leav- 
ing the Luttinger parameters as adjustable parameters, 
or only one correlation function could be measured in a 
given experiment, not allowing to test for universality of 
the description. 

We also gave a detailed description of the dynamical 
spin-spin correlations, for T = using the time depen- 
dent DMRG method, for a wide range of energies and all 
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FIG. 22. (Color online): Inelastic neutron scattering inten- 
sity measured along a* of BPCB^ with a tt momentum in 
the rungs (Q • di = tt) at T = 250 mK after substraction of 
the zero- field background. In each panel, cuts at fixed en- 
ergy (shown by white dashed lines in Fig. [20]) are plotted: 
(a) uj - 0.4 meV and rn z = 0.24, (b) uj - 0.4 meV and 
m z = 0.72, (c) uj = 0.2 meV and m z = 0.24, (d) uj = 0.2 meV 
and m z = 0.72. The circles correspond to the experimental 
data. The solid red (black) curves are the theoretical data 
for the ladder (the spin chain mapping) convolved with the 
instrumental resolution. The shaded bands indicate the error 
bar in the experimental determination of a single proportion- 
ality constant valid for all fields, energies, and wave vectors. 
The width of these areas combines the statistics of all our 
scans with uncertainties in the exact magnetization values at 
the chosen fields and in the convolution procedure. 



momenta. The excitations reveal many important infor- 
mations on the system and are well suited to characterize 
it. In particular we show the interesting evolution of the 
excitations in the system with the magnetic field and with 
different coupling strengths. Quite interestingly the in- 
termediate energy part can be related to the excitations 
of a t- J model and shows thereby features of itinerant sys- 
tems. We also showed that the dynamical correlations of 
the ladder posses characteristic high energy features that 
are clearly distinct from the corresponding spectrum for 
spin chains. 

The numerical calculation is efficient for the high and 
intermediate energy part of the spectrum for which the 
Luttinger liquid description cannot be applied. We 
showed that the two methods, numerics and LL have 
enough overlap, given the accuracy of our calculation so 
that we can have a full description of the dynamical prop- 
erties at all energies. This allowed us to use each of the 
method in the regime where it is efficient. In particu- 
lar, in the present study we did not push the numerical 
calculations to try to obtain the exact behavior at low 



energies, but focus on the high and intermediate energy 
regime. We used the analytical description coupled to the 
numerical determination of the Luttinger parameters to 
obtain a very accurate low energy description. We made 
the connection between our results and several analyti- 
cal predictions. In particular at intermediate- low energy 
our calculation agrees with the Luttinger liquid predic- 
tion of incommensurate points and behavior (divergence, 
convergence) of the correlations. 



We compared our numerical results with existing INS 
data on the compound BPCB and found excellent agree- 
ment. It is rewarding to note that the resolution of our 
dynamical calculation is, in energy and momentum, at 
the moment better than the one of the experiment. The 
comparison between theory and experiment is thus essen- 
tially free of numerical error bars. Given the current res- 
olution of the INS experiment it is difficult to distinguish 
in the low energy part of the spectrum the difference be- 
tween the dynamical correlation of the true ladder and 
the one of an anisotropic spin 1/2 system, which corre- 
sponds to the strong rung exchange limit. More accurate 
experiments would be desirable in this respect. An alter- 
native route is to probe experimentally the high energy 
part of the spectrum, since these high energy excitations 
contain many characteristic features of the underlying 
model. 



On the conceptual side and also in connection with 
the compound BPCB several points remain to be inves- 
tigated. An extension of the dynamical results to finite 
temperature would be desirable. This could be used to 
study different effects such as the interesting shifts and 
damping of the triplet modes with temperature that have 
been observed in three dimensionnal gapped system^. A 
second important step would be to improve the descrip- 
tion of the quasi one dimensional systems, by including 
in a mean-field way the effect of the other ladders in the 
numerical study. This is something we already partly 
performed, but the extension to the dynamical quantities 
remains to be done. This is specially important close to 
the quantum critical points h c \ and h C 2 where the inter- 
ladder coupling becomes crucial and the system under- 
goes a dimensional crossover between a one dimensional 
and a higher dimensional (three dimensional typically) 
behavior. Understanding such a crossover is a particu- 
larly challenging question since the system goes from a 
description for which an image of essentially free fermions 
applies (in the one dimensional regime) to one for which 
a description in terms of essentially free bosons (the three 
dimensional regime) applies. 



One step further would be the extension of the investi- 
gation of the combined numerical-analytical methods we 
used here to other systems including ladder systems with 
certain asymmetries or in the presence of disorder. 
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Appendix A: Strong coupling expansion of a single 
spin ladder 

In this appendix we show how the spin ladder Hamil- 
tonian (J2j) at strong coupling (7 < 1) can be expressed 
in bosonic operators acting on single bonds introduced 
in Ref. [94]. This representation allows a classification of 
the excitations due to their position in energy. We first 
derive perturbatively an effective system based on this 
Hilbert space organization by energy sectors. We intro- 
duce the Schrieffer- Wolff transformation that maps the 
physical system to the effective, and approximate it us- 
ing a strong coupling expansion. Thus, we evaluate the 
rung densities of the ground state in the spin liquid, and 
derive an effective theory for the gapless regime. Fur- 
thermore we evaluate the deviation of the LL parameters 
from the spin chain mapping. 



The Hermitian operator A can be expanded in powers of 

7 

A = Ai + 7A2 + • • • . (A4) 
Thus H e ff can be written in orders of 7 as 

J^Hes = H ± + + j 2 HW + . . . , (A5) 

where 

HW=H\\+i[A 1 ,H ± ], (A6) 
HW=i[A 2 ,H ± ] - ±[A U [A U H ± ]] + i[A u H\\], (A7) 

etc. Through this expansion, the unitary transform e LlA 
can be perturbatively determined computing the A k re- 
cursively in order to eliminate the transitions between 
the energy sectors of excitations in H e ^. Since the first 
term J±H± in Eq. (|A5p leads to a separation of excita- 
tions on the order of the energy scale J± (cf. Fig.(2ja) the 
decoupled bond limit provides the effective Hilbert space 
that contributes to each energy sector. The second term 
J\\H^ causes broadening of these bands on the order of 
J|| and can induce a complex structure within the energy 
bands. To obtain the desired expansion up to the first 
order in 7 we choose 

Ai = j y^44+i (^,0^+1,0 - u,+u+i- 
1 

-tz,-*z+i,+)+h.c, (A8) 
where h.c. stands for the Hermitian conjugation. 



1. Strong coupling expansion 

The four-dimensional Hilbert space on each rung I is 
spanned by the states and \t~), (cf. Eqs. (j9j) 

and ([TO]) ), obtained by applying the boson creation oper- 
ators sj, + , t\ , and t\ _ to a vacuum state. A hardcore 
boson constraint applies on each rung /, i.e. 

Ql,s + Ql,+ + Qi,o + Ql- = 1 (AI) 

where qi j3 = sjsi and g^ k = 4/A&> k = ±,0. 

While the Hamiltonian on the rung H± ((3]) is quadratic 
in the boson operators 

L 3 
H± = X] [(1 - h z /J±)gi,+ + ao + (1 + h z /J ± ) ei ,-]-~L, 
1=1 

(A2) 

the chain Hamiltonian H\\ dU) is quartic, and its structure 
is quite complex. The advantage of the boson represen- 
tation reveals when considering the case of small 7. In 
that case we perform a Schrieffer- Wolff transformation of 
the spin ladder Hamiltonian (|2j) 

H eE = e^ A He~^ A . (A3) 



2. Rung state density in the spin liquid 

In the spin liquid phase, the decoupled bond limit pro- 
vides the effective ground state |O e ff) = \s • • • s) which is 
related to the physical ground state by 

|0)=e- < ^|0 eff ). (A9) 

So the triplet density of the ground state (pk) (with k = 
=t, 0) is given by 

(Pk) = (m,k) = (*••• s\e^ A Quk e-^ A \s ■■■s). (A10) 

Using Eq. (|A8|) , and keeping only the non- vanishing cor- 
responding terms in (|A10|) up to second order we get 

{pk)=l 2 {s---s\A lQl , k A 1 \s---s) = (All) 

o 

In the case of the compound BPCB (see Sec. IVI A[) this 
expansion gives {pk) = 0.01, and due to the hardcore 
boson constraint (Eq. (|Al]) ) (p s ) = (g^ s ) = 0.97. Even 
though we took into account only the first order term for 
A in Eq. (|A4[) . this approximation of the triplet density 
differ from the direct numerical computations (in Fig. [4]) 
of only - 20%. 
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3. Effective Hamiltonian in the gapless regime 

The first order term of the effective Hamilto- 

nian (|A5|) is computed substituting (|A8|) into (|A6|) . This 
leads to in the form 



fc=0 z 



(i) 



(A12) 



where # ± is given by Eq. (lAT9|) and by Eq. (fAT3|) . 
Following Ref. Q3, we map the two low energy modes on 
the two states of a pseudo spin (Eq. (fTTj) ) and replace the 
boson operators and t+ with the spin- 1/2 operators 



S, + = 5i-=a,V, 5f = e,, + --. (A21) 



The effective Hamiltonian is the spin- 1/2 XXZ Heisen- 
berg Hamiltonian, Eq. ([T3]) . 



where 



^o!z = + 7^z+i,+£z,+ + h - c -, ( A1 3) 



= 5 I+i4,o^+i,o^z+^ T + i,+^ T ,o^+i,o^,+ + h-c, (A14) 



t *t 



+■ 4+i,o4,o^+i,-^,+ + 4+i,o4,o^+i,+^,- + hx -> ( A15 ) 



- 4+i,o4,-^+i,-^,o + h.c, 



and 



ft 



(i) 
4,/ 



Ql+i-Ql,- 



(A16) 



(A17) 



Here we regrouped the terms such that each J±H^ acts 
on the corresponding energy sector k J_l, k = 0, 1, . . . , 4 in 
the gapless regime. Note that in each sector Ai = such 
that to the considered order in 7 the Hamiltonian (j2j) 
corresponds to the effective Hamiltonian. 



b. Sector of energy J± 

The effective Hilbert space of the J± energy sector cor- 
responds to a single \t°) triplet excitation lying in a sea of 
singlets \s) and triplets The effective Hamiltonian 

up to first order in 7 is given by 



H = J^H ± + J,, (h™ + H[ V ' 



(A22) 



The excitation \t°) can be interpreted as a single hole 
excitation in a spin chain formed by \s) and Each 
rung state of this sector is identified with 



li> = l*>, lt) = |t + ), |0) = |t°). 



(A23) 



In this picture the Hamiltonian ()A22|) can be mapped 
onto the anisotropic t-J model 



Ht-J — H: 



xxz 



H t + H s . h + e. 



where e = (J± + h z )/2 is an energy shift. 
The hopping term 



C 3"+l,cr C Z,(J 



(A24) 



(A25) 



a. Low energy sector 

When focusing on the low energy sector, the \s) and 
the \t+) modes dominate the behavior and we can assume 
a vanishing density of \to) and triplets. Thus the 
hardcore boson constraint (|A1|) simplifies to 



(A18) 



Qi,8 + Qi,+ = 1 
and the rung Hamiltonian (|A2|) to 



Hx = (l-^x)E«.+ -f i - (A19) 

Further the only contribution to the first order term in 
7 comes from Hq 1 ^. Taking this into account we obtain 
from Eq. (|A5|) for the Hamiltonian ((2]) in the lowest en- 
ergy sector 



H=J ± H ± + J ll H^ 1 \ 



stems from the term J\\H[ in Eq. (jA22|) . Here cj a (q ?0 -) 
is the creation (annihilation) operator of a fermion with 
pseudo spin a = f, ^ at the site I. Note that although 
we are dealing here with spin states, it is possible to 
faithfully represent the three states of each site's Hilbert 
space (|s), \t )) using a fermion representation. 

Additionally a nearest neighbor density-density term 
between the up spin and the hole arises 



(A26) 



(A20) 



Here n^h is the density operator of the hole at site I. This 
term stems from the nearest-neighbour interaction be- 
tween the triplets, i.e. the second term in Equ. (|A13|) . 
Mapping this term onto a spin chain in the presence of 
a hole leads to an interaction term between the hole and 
the spin-up state. Note that this is in contrast to the 
usual mapping on a spin chain without a hole in which 
case it would only causes a shift in energy. 
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4. Second order perturbation and Luttinger liquid 
parameters 

The second order term H^ (|A7|) of the expan- 
sion (|A5|h contains a huge amount of terms. Neverthe- 
less considering the low energy sector, only the following 
terms 

„.(2) 3 
= ~o / ,Ql,sQl+i,s 

1 

\ + h - c -) ( A27 ) 



are important. The first term in Eq. (jA27|) is a singlet 
density-density interaction which can be absorbed into 
the coupling of the XXZ chain, and the second term is a 
conditionnal hopping^. In order to study the effects of 
Hq 2 ^ on the LL parameters u and K (see Fig. [23]), we first 
replace the boson operators with the spin- 1/2 operators 
(Eq. (|A21)0 . So the Hamiltonian ((2]) in the low energy 
sector becomes 



H — Hxxz - 



1 



J i-i 



^1+1 



h.c. 



+const 



(A28) 



where Hxxz is the XXZ Heisenberg chain Hamiltonian 
(Eq. (fl3j) ) with the corrected parameters 



A (2) 



77 



h 



r,( 2 ) 



h z 



(A29) 



up to the second order in 7. For the BPCB parame- 
ters (Eq. d5j) A( 2 ) ^ 0.4 instead of A = 0.5 for the spin 
chain mapping (first order approximation). The LL pa- 
rameters u, K and A x of Hxxz with the anisotropy A^ 2 ^ 
are computed, and we treat the three terms interaction 
(conditional hopping) approximating 1/2 — Sf = 1/2— m z 
(mean- field approximation). The remaining term is then 
bosonized using the expression ([TT)) for S ± (x = I). It 
leads to the corrected LL parameter u and K of the 
Hamiltonian (|A28|) through the relations 



the low energy physics of the gapless regime in the spin 
ladder model ([2]) and the spin chain mapping (fT3|) . We 
start with the determination of the parameters that to- 
tally characterize the LL description. Then we summa- 
rize the spin-spin correlation functions deduced from the 
LL theory at zero and finite temperature. 



1. Luttinger liquid parameter determination 

In this appendix we detail the determination of the 
LL parameters u, K and the prefactors A x , B x and A z 
(see Eqs. (HU), (jT7|) and (fT5|) V The parameters K, A x , 
B x and A z and their dependence on the magnetic field 
have been previously determined in Refs. l25[|26l andl96l 
for different values of the couplings than considered here. 
We obtain these parameters in two step o 3Q i 97 : 

(i) We determine the ratio u/K from its relation to the 
static susceptibility, Eq. (|B1|) , which is a quantity 
numerically easily accessible with DMRG or Bethe 
ansatz (for the spin chain only). 

(ii) The parameters K and the prefactors A x , B x and 
A z are extracted by fitting numerical results for 
the static correlation functions obtained by DMRG 
with their analytical LL expression. For the spin 
chain, it is also possible to deduce the product 
uK from the magnetic stiffness computed by Bethe 
ansatz&2£. 

We give the relations used for the spin chains only, since 
from these the relations for the spin ladders can be easily 
inferred using the following relations 



uK = uK + 27r~f 2 J ± A X (1/2 -rh z ) 
u/K = u/K 



(A30) 



The corrected u and K are plotted in Fig. [231 and clearly 
show the asymmetric signature of the full ladder param- 
eters induced by the conditional hopping term in (jA28p . 
Note that the lack of convergence K 1 when m z = 
m z + 1 1 2 —> is obviously an artefact of the mean- field 
approximation 1/2 — Sf = 1/2 — m z and the big sensitiv- 
ity in the K determination with the ratio of Eq. (|A30|) . 



dm z 



drh z 
dh z 



(Sf 9V Sf t9V )^2(SfSf,) , <S* ) <5f ) + - 
(Sf,oSf, >) -+ {SfSf) + \ ((Sf) + (Sf,)) + \. 



a. Susceptibility and u / K 



The LL theory predicts that the static susceptibility 
is related to the ratio u/K through the relation^^ 

dh z 1 



U 

K 



^ drh z 



(Bl) 



We numerically compute the static susceptibility using 
DMRG or Bethe ansatz (for the spin chain only) and 
infer the ratio u/K with a negligible error. 



Appendix B: Luttinger liquid 



b. Static correlation functions 



In the following, we present several properties of the 
Luttinger liquid providing a quantitative description of 



For the extraction of the parameters K, A x , B x and 
A z , we fit numerical DMRG results for the correlation 
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FIG. 23. (Color online): LL parameters u, K and the prefac- 
tors of the spin operators A x , B x , A z versus the magnetiza- 
tion per rung m z computed for a spin ladder with the BPCB 
couplings (|39p (red crosses) and for the spin chain mapping 
(blue stars) . The strong coupling expansion of u and K up to 
the second order in 7 (discussed in appendix I A 4|) is plotted 
in black dashed lines. 



functions (SfSf,), (SfSf) and the local magnetization 
(Sf) with their analytical expression for finite system 
size, Eqs. (11), (12) and (13) in Ref. I26I respectively. In 
the limit of infinite system size and far from the bound- 
aries these LL correlations simplify to the well known 
power law decay 



we use the previously extracted value for K to fit the 
longitudinal correlation (zz-correlation (SfSf)) and the 
magnetization, (Sf), which allow us to determine A z . 
The values determined by both fits are very close and in 
Fig. [23] the average value of both is shown. 

All the results presented in Fig. [23] were done for L = 
200 and several hundred DMRG states after an average 
on the four sets of used data points in the fit 10 < Z, V < 
170, 30 < 1,1' < 170, 10 < 1,1' < 190, 30 < 1,1' < 190. 
The error bars correspond to the maximum discrepancy 
of these four fits from the average. We further checked 
that different system lengths lead to similar results. 

The LL parameters of the ladder system ((2]) for the 
BPCB couplings ([39]) are presented in Fig. [23] as a func- 
tion of the magnetization per rung. Additionally we show 
the parameters of the spin chain mapping (computed for 
the spin chain Hamiltonian ([13]) ) for comparison. When 
the ladder is just getting magnetized, or when the ladder 
is almost fully polarized, K 1 (free fermion limit) and 
u —> (because of the low density of triplons in the first 
case, and low density of singlets in the second case). For 
the spin chain mapping, the reflection symmetry around 
m z = 0.5 arises from the symmetry under ir rotation 
around the x or y— axis of the spin chain. This sym- 
metry has no reason to be present in the original ladder 
model, and is an artefact of the strong coupling limit, 
when truncated to the lowest order term as shown in ap- 
pendix [A] The values for the spin ladder with the com- 
pound BPCB parameters can deviate strongly from this 
symmetry. The velocity u and the prefactor B x remain 
very close to the values for the spin chain mapping. In 
contrast, the prefactors A z , A x and the exponent K de- 
viate considerably and A x and K become strongly asym- 
metric. The origin of the asymmetry lies in the contri- 
bution of the higher triplet states 2 -, and can be under- 
stood using a strong coupling expansion of the Hamilto- 
nian ([2]) up to the second order in 7 (see appendix IA 4j) . 
This asymmetry has consequences for many experimen- 
tally relevant quantities and it was found to cause for 
example strong asymmetries in the 3D order parameter, 
its transition temperature and the NMR relaxation rate 
(see Fig. [TUJ Fig. M and Fig. [8]). 



(SfSf,) = A x 



(-1) 



l-V 



\l - 



B x (~l) 



l-V 



cos[q(l - Q] 



(B2) 



<5f5?>=m« a + ^(-l)'-'' 00B[9(/ - r)] 



K 



\l-l'\ 



i\1K 



2ir 2 \l-l'\ 2 
_(B3) 

and the local magnetization becomes constant, (Sf) = 

We first fit the transverse correlation (^-correlation 
(SfSfr)) to extract the parameters A x , and B x . Then 



2. Dynamical Luttinger liquid correlations 



In this appendix we summarize previous results of the 
Luttinger liquid description of the dynamical correla- 
tion function o 19 i 2Q i" at zero temperature. Note that the 
weak and strong coupling limits which have been studied 
separatel y 19 ! 20 can be connected 2 . The expression of the 
correlations in the whole regime is given by^ 
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k 2 A z 
uT(K) 2 



+ e(u - u\q - 2tt(1 -m z )\) 



Q(uj — u\q — 2i\m z \ 

4^ 

uj 2 -u 2 (q- 27r(l-ra 2 )) 2 



4u 2 



uj 2 — u 2 (q — 2irm z ) 2 

\-K' 



1-K 



Kuj 

+ — ©M ~ uq) + S(oj + ug)] (B4) 



Stt 2 A x 



uT(l/4K) 2 



6(cj — — 7r| 



cj 2 — u 2 (q — 7r) 5 



1-1/4K 



+ 



4tt 2 B x 



0(uj - u\q - tt(1 -2m z )\) 



uT(rj + )T(rj_) 

+G(u-u\q- 7r(l + 2m*)|) 



2ix 



fc[g-7r(l -2m z )] 



l-rj- 



2u 



i[q - tt(1 + 2m 2 



1-77+ 



2u 



uj + u[g — 7r(l — 2m 2 )] 



1-77+ 



2ix 



cj + ?i[g-7r(l + 2ra 2 )] 



1-77- 



(B5) 



* (a) 




tt(1 - 2m 2 ) 7T tt(1 + 2m' 




tt(1 - 2m 2 ) 7T tt(1 + 2m 2 ) 



FIG. 24. (Color online): Map of the momentum- frequency 
low energy correlations of the LL model where the white ar- 
eas represent the continuum of excitations. In the striped 
areas no excitations are possible, (a) S§ z (q,uj): the dash- 
dotted lines (blue) are the excitation peaks close to q = 0, 2ty 
and the dashed lines (red) are the continuum lower boundary 
with edge exponent 1 — K close to q — 27rm z ,27r(l — m z ). 
(b) Sn~(q,w), (c) S~ + (q,uj): the continuum lower boundary 
close to q = 7r, 7r(l =b 2m z ) is represented by solid lines (black) 
(edge exponent 1 — 1/4K), dashed lines (red) (edge exponent 
l—ry- — 2 — 1/4K — K) and dash-dotted lines (blue) (edge 
exponent 1 — 77+ = — 1/ 4K — K) . 



with the Heaviside function Q(x) and rj± = 1/4K±1+K. 
The correlation S~^~ is obtained replacing m z — >• — m z in 
the S+~ expression Eq. (|B5|) . 

The expressions Eq. (|B4[) and Eq. (|B5[) exhibit the 
typical behavior of the frequency- momentum LL correla- 
tions: a continuum of low energy excitations exists with 
a linear dispersion with a slope given by the Luttinger 
velocity ±iz. The spectral weight at the lower boundary 
of the continuum displays an algebraic singularity with 
the exponents related to the Luttinger parameter K. A 
summary of this behavior is sketched in Fig. [2U For the 
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FIG. 25. (Color online): (a) Different exponents that appear 
in the LL correlation functions, Eqs. (|B4|) and (|B5|) . versus the 
magnetization m z . The solid (dashed) lines are determined 
from the ladder (spin chain mapping) exponent K in Fig. [23l 
The exponent 1 — K of the S zz correlations is shown in (a.l), 
and the exponent 1 — 1/4K of the correlations at the 

q — 7r branch in (a. 2) (lower red curves). The exponents 
1 — 77- = 2 — 1/4K — K (upper black curves) in (a.2) and 
1 — 77+ = — 1/4K — K in (a.3) correspond of both sides of the 
incommensurate branches of the (see Fig. [24]) . (b) Shift 
of the ground state energy per rung versus the magnetization 
AE (m z ) = E (0) - E (m z ). 



considered system the longitudinal correlation S§ z is pre- 
dicted to diverge with the exponent 1 — K at its lower 
edge. As shown in Fig.[25]a.l the exponent of this diver- 
gence is very weak < 0.2 for the parameters of BCPB. 
The transverse correlations exhibit a distinct be- 

havior depending on the considered soft mode. Close 
to q = 7r the weight diverges with an exponent given by 
1 — 1/ 4K. This divergence is strong for the considered pa- 
rameters {1-1/ w 3/4 > in Fig.[25ja.2). In contrast 



31 



at the soft mode q = 7r(l — 2m z \ 7r(l J r 2m z ) a divergence 3. Finite temperature transverse staggered 

(cusp) is predicted at the lower edge with the exponent Luttinger liquid correlation 

2 - 1/4K - K w 3/4 in Fig. E3a.2 (-1/4X — -5/4 

in Fig. [25la.3). The determination of the relaxation time T-f 1 (see 

Sec. IIV C[) and of the transition temperature T c to the 
3D-ordered phase (see Sec. lIVDT]) requires the staggered 
transverse retarded correlation function 

X x a x (x, t) = -ie(t)(-lf ( [S x (x, t), S x (0, 0)] ) (B6) 

written in term of the spin chain mapping operators (fl2j) 
with 6(t) the Heaviside function. In the gapless regime, 
using the bosonization formalism ([I7|) . and taking into 
account only the most relevant term, we can compute it 
as described in Ref. [6| for the LL Hamiltonian (IT6|) : 



X X a X (x,t) 



-e(t)Q(ut - x)0(ut + x) — 



2A x sm{^) 



(B7) 



and its Fourrier transform: 

Ac sin (^) /2tt\^ 



(cj — wg) 

47T 



47T 



J-i-J- 

81T' 4.K" . 



(B8) 



where £?(£, ?/) = ^ 



Appendix C: DMRG method 

In this appendix, we describe first the time dependent 
DMRG method and its extension to finite temperature. 
Then we give the technical details related to the compu- 
tation of the momentum- frequency correlations. 



1. Time dependent DMRG 

The t-DMRG^Zr— (time dependent DMRG) method is 
based on the principle of the original DMRG to choose 
an effective reduced Hilbert space to describe the physics 
one is interested in. Instead of choosing only once the 
effective description for the evolution of the state, the 
t-DMRG adapts its effective description at each time- 
step. The implementation of this idea can be performed 
using different time-evolution algorithms. Here we use 
the second order Trotter- Suzuki expansion for the time- 
evolution operator using a rung as a uni t 38 i 39 . The errors 
arising in this method are the so-called truncation error 
and the Trotter - Suz uki error. Both are very well con- 
trolled (see Ref. llOOl for a detailed discussion). 



2. Finite temperature DMRG 

The main idea of the finite temperature DMRG^— 
(T-DMRG) is to represent the density matrix of the phys- 
ical state as a pure state in an artificially enlarged Hilbert 
space. The auxiliary system is constructed by simply 
doubling the physical system. Starting from the infinite 
temperature limit the finite temperature is reached evolv- 
ing down in imaginary time. The infinite temperature 
state in this auxiliary system corresponds to the totally 
mixed Bell state |^(0)) = H1L1 Ea z where 

I Xi Xi) is the state at the bond I of the auxiliary system 
that has the same value Xi on the two sites of the bond 
(physical and its copy). The sum is done on all these 
N\ states |A/). We evolve the physical part of |^(0)) i n 
imaginary time to obtain 

\m) = e- pH/2 \m) (ci) 

using the t-DMRG algorithm presented in appendix IC II 
with imaginary time. We renormalize this state at each 
step of the imaginary time evolution. Thus the expecta- 
tion value of an operator O acting in the physical system 
with respect to the normalized state is directly 

related to its thermodynamic average, i.e. 

{o) P = = (m\o\m) (C2) 

at the temperature T = We use this method to 

compute the average value of the local rung magnetiza- 
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tion m z (T) and energy per rung E(T) in the center of the 
system. Additionally we extract the specific heat c(T) by 

cC8 + 5(3/2) « - (/? + ^ /2)2 {(E)p + $p - (E)p) (C3) 

where (5/3 is the imaginary time-step used in the T- 
DMRG. 

To reach very low temperatures T for the specific 
heat, we approximate the energy by its expansion in T 

n 

E(T)^E Q + Y. a * Tl ( C4 ) 

i=2 

up to n = 4. The energy at zero temperature Eq is 
determined by a standard T = DMRG method. Since 
E(T) has a minimum at T = the linear term in the 
expansion (|C4p does not exist. The numbers oli (i = 
2,3,4) are obtained fitting the expansion on the low T 
values of the numerically computed E(T). 

Typical system lengths used for the finite temperature 
calculations are L = 80 (L = 100 for the spin chain map- 
ping) keeping a few hundred DMRG states and choosing 
a temperature step of 6/3 = 0.02 K" 1 (6/3 = 0.01 K" 1 for 
the spin chain mapping). 

Let us note that recently a new method has been de- 
veloped to treat finite temperature o 101 i 102 which is very 
promising to reach even lower temperatures. 



3. Momentum-frequency correlation functions 

To obtain the desired frequency-momentum spin-spin 
correlations (SqP(q,u) in (|3T|) ). we first compute the cor- 
relations in space and time 

Stfitn) = (0|^^^ + L/2,^-^^^/ 2 , 1 |0) (C5) 

with I = -L/2 + l,-L/2 + 2, ...,L/2, k = 1,2, and t n = 
n6t (n = 0, 1, . . . , N t ) is the discrete time used. These 
correlations are calculated by time-evolving the ground 
state |^o) = |0) and the excited state = S^/ 21 \0) 
using the t-DMRG (see appendix IC 1|) . 

Afterwards the overlap of S^_ L , 2 k \ipi (t)) and \^o(t)) 
is evaluated to obtain the correlation function (|C5|) . 

In an infinite system reflection symmetry would be ful- 
filled. To minimize the finite system corrections and to 
recover the reflection symmetry of the correlations, we 
average them 

\ (S a _{ k {t n ) + Sf k (t n )) Sf k (t n ). (C6) 

We then compute the symmetric (antisymmetric) corre- 
lations (upon leg exchange) (see Sec. IV A[) 

Sf qy {t n ) = 2(S%(t n ) ± S$(t n )) (C7) 



with the rung momentum q y = 0, 7r respectively. Finally, 
we perform a numerical Fourrier transform^ 

S^(q,w)»St f; Yl e^-^Sf qy (t n ) 

n=-N t + l l=-L/2+l 

(C8) 

for discrete momenta q = 2irk/L (k = 0, 1, . . . , L — 1) and 
frequencies Co 104 . The momentum q has the reciprocal 
units of the interrung spacing a (a = 1 is used if not 
mentioned otherwise). Due to the finite time step 6t, 
our computed SqP(q,uj) are limited to the frequencies 
from —7r/6t to ir/6t. The finite calculation time tf = 
N t 6t induces artificial oscillations of frequency 2it/tf in 
SqP(q,uj). To eliminate these artefacts and reduce the 
effects of the finite system length, we apply a filter to 
the time-space correlations before the numerical Fourrier 
transform (|C8)h i.e. 

S?£(t n )f{l,t n )^S?l(t n ). (C9) 

We tried different functional forms for the filter f(l,t n ) 
(cf. Ref. [38| as well). In the following the results are ob- 
tained by a Gaussian filter f(l,t n ) = e -(^/L) 2 -{2t n /t f f 
if not stated otherwise. As the effect of this fil- 
tering on the momentum-energy correlations consists 
to convolve them by a Gaussian function f(q,uj) = 
t/L/(327r)e" (u; ^ /4)2 " (<?L/8)2 , it minimizes the numerical 
artefacts but further reduces the momentum-frequency 
resolution. 

Typical values we used in the simulations are system 
lengths of up to L = 160 sites while keeping a few hun- 
dreds DMRG states. We limited the final time tf to be 
smaller than the time necessary for the excitations to 
reach the boundaries (tf ~ L/2u with u the LL velocity 
in Fig. [23]) in order to minimize the boundary effects. The 
computations for the BPCB parameters, Eq. ([39]) . were 
typically done with a time step of 6t = 0.0355 J\\~ x up to 
tf = 71 J\\~ x (but calculating the correlations only every 
second time step). The momentum- frequency limitations 
are then 6uo w 0.11 Jn and 6q ~ 0.1 a -1 . Concerning the 
other couplings and the spin chain calculations, we used 
a time step St = 0.1 Jy -1 up to t/ = 100 Jy -1 (also 
with the correlation evaluations every second time steps) 
for a momentum- frequency precision 6uo w 0.08 Jn and 
6q w 0.1 a~\ 

Different techniques of extrapolation in time (using lin- 
ear prediction or fitting the lon g tim e ev oluti on with a 
guess asymptotic form cf. Refs. TToBI and Il06l ) were re- 
cently used to improve the frequency resolution of the 
computed correlations. Nevertheless, as none of them 
can be apply systematically for our ladder system due to 
the presence of the high energy triplet excitations (which 
result in very high frequency oscillations), we decided not 
to use them. 
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Appendix D: Quantum Monte Carlo determination 
of the 3D-ordering transition T c 

The three dimensional network of couplings of the cou- 
pled ladder Hamiltonian Eq.[T]and shown in Fig.[TH]is not 
frustrated, and can therefore be simulated using Quan- 
tum Monte Carlo algorithms. We employ a stochastic 
series expansion implementation with directed loop up- 
dates 64 provided with the ALPS libraxies25*2£. This nu- 
merical exact approach is complementary to mean-field 
approaches, because the later tend not to be quantita- 
tively accurate due to the neglect of quantum fluctuations 
in the interladder coupling. 

In order to detect the transition temperature we mea- 
sure the spin stiffness p s based on winding numbe rs in 
the three spatial directions. As pointed out in Ref. llQTt 
when plotting Lp s for different system sizes L, the dif- 
ferent curves cross at T c . Alternatively we measured the 
squared order parameter w? x . The quantity L 1+7? ra^ (as- 
suming the 3D XY universality class value of r] ~ 0.04) 
also crosses at T c when plotted for different L. 

Previously the 3D ordering temperatures of coupled 
spin ladders in a magnetic field have been determined 
using a specific feature of the magnetization m(T)2£. It 
turns out, that while the feature in m(T) indeed locates 
T c accurately for simple coupled dimer systems^, the 
same does not hold for coupled ladders. In the ladder 
case one has to resort to the spin stiffness crossing or 



order parameter measurements to locate T c accurately. 

When simulating the coupled ladder Hamiltonian, a 
suitable finite size sample setup is required. Due to the 
spatial anisotropics in the problem - the ladder direction 
being singled out over the two transverse spatial direc- 
tions - an appropriate aspect ratio needs to be kept 72 . 
In our simulations we chose an aspect ratio of about 12, 
i.e. the samples were 12 times longer along the ladder 
direction than in the transverse directions. 

In Fig. [26] we show simulation results for one partic- 
ular set of couplings: the rung and leg couplings were 
set to J± = 12.9 K and J\\ = 3.3 K respectively, the 
^-factor was 2.17, the magnetic field amounted to 8.9 T 
and the interladder coupling J' was set to 80 mK. These 
couplings are identical to those used in Ref. [3l|. In the 
upper panel we show the rescaled spin stiffness Lp s for 
two different system sizes (768 versus 6144 spins). One 
locates a crossing at 210 mK for this observable. In the 
lower panel we show the rescaled squared order parame- 
ter L 1+7? m 2 , which also exhibits a crossing at the same 
temperature. These matching estimates for the critical 
temperature render us confident that we accurately lo- 
cate the critical temperature. 

Based on this and subsequent simulations either with 
an identical J' but a higher magnetic field of 11.9 T, 
or the same field and a smaller J' = 50 mK, we are 
able to determine and verify the use of single rescaling 
factor a = 0.74(1) relating the real and effective mean- 
field interladder coupling^ 2 -, as presented in Sec. IIVD 1[ 
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